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Abstract 


A  control  chart  is  often  used  to  detect  a  change  in  a  process.  Following  a  control 
chart  signal,  knowledge  of  the  time  and  magnitude  of  the  change  would  simplify  the 
search  for  and  identification  of  the  assignable  cause.  In  this  research,  emphasis  is  placed 
on  count  processes  where  overdispersion  has  occurred.  Overdispersion  is  common  in 
practice  and  occurs  when  the  observed  variance  is  larger  than  the  theoretical  variance  of 
the  assumed  model.  Although  the  Poisson  model  is  often  used  to  model  count  data,  the 
two-parameter  gamma-Poisson  mixture  parameterization  of  the  negative  binomial 
distribution  is  often  a  more  adequate  model  for  overdispersed  count  data.  In  this  research 
effort,  maximum  likelihood  estimators  for  the  time  of  a  step  change  in  each  of  the 
parameters  of  the  gamma-Poisson  mixture  model  are  derived.  Monte  Carlo  simulation  is 
used  to  evaluate  the  root  mean  square  error  performance  of  these  estimators  to  determine 
their  utility  in  estimating  the  change  point,  following  a  control  chart  signal.  Results 
show  that  the  estimators  provide  process  engineers  with  accurate  and  useful  estimates  for 
the  time  of  step  change.  In  addition,  an  approach  for  estimating  a  confidence  set  for  the 
process  change  point  will  be  presented. 
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CHANGE-POINT  METHODS  FOR  OVERDISPERSED  COUNT 

DATA 

1.  Introduction 

This  chapter  will  start  with  some  relevant  background  information  on  statistical 
process  control  (SPC)  that  sets  the  foundation  for  this  research  effort.  Following  the 
background  review,  the  problem  definition  will  be  introduced  along  with  the  research 
objectives  and  assumptions  of  this  thesis.  The  last  section  of  this  chapter  will  lay  out  the 
organization  of  this  thesis. 

1.1  Review  of  Statistical  Process  Control 

The  focus  of  this  research  is  on  change-point  estimation  within  an  SPC  context. 
Before  going  further  with  this  research,  a  brief  overview  of  SPC  topics  and  terminology 
pertaining  to  this  research  is  reviewed.  Refer  to  Montgomery  [9]  and  Ryan  [19]  for  more 
detailed  descriptions  of  SPC  and  its  applications. 

1.1.1  SPC  Terminology  and  Concepts 

Montgomery  [9]  defines  statistical  process  control  as  a  set  of  problem-solving 
tools  that  can  be  applied  to  a  process  to  achieve  stability  and  improve  capability  by 
reducing  variability.  This  tenn,  process,  is  usually  thought  of  in  a  manufacturing  setting, 
but  can  include  non-manufacturing  settings  as  well.  Non-manufacturing  processes  can 
include  tracking  errors  in  some  paperwork  process  of  a  business,  or  even  monitoring 
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daily  coalition  casualties  in  the  Iraq  war,  which  is  examined  in  Chapter  4.  Of  the 
different  SPC  tools  available,  such  as  histograms  and  cause-and-effect  diagrams,  this 
research  will  only  focus  on  the  control  chart,  specifically  the  Shewhart  control  chart  and 
the  cumulative  sum  control  chart.  Before  getting  into  the  specifics  of  each  of  these 
control  charts,  some  overall  SPC  terminology  needs  to  be  addressed. 

Any  process,  no  matter  how  finely  tuned,  will  always  contain  some  variability. 
Natural  variability,  or  uncontrollable  variability,  is  inherent  in  any  process.  A  process 
operating  in  a  state  only  affected  by  natural  variability  is  said  to  be  statistically  in¬ 
control.  However,  other  kinds  of  variability  may  be  present  in  the  output  of  a  process. 
For  instance,  in  a  manufacturing  context  these  non-natural  sources  of  variability,  or 
controllable  sources  of  variability,  may  be  introduced  into  a  process  by  improperly 
calibrated  machines,  defective  raw  materials,  or  even  human  error.  Causes  of  these 
sources  of  variability  are  referred  to  as  assignable  causes.  A  process  is  said  to  be  out-of¬ 
control  if  assignable  causes  are  present.  Thus,  one  of  the  primary  objectives  of  SPC  is  to 
quickly  detect  a  change  in  the  process,  look  for  assignable  causes,  and  take  corrective 
action  to  bring  the  process  back  to  a  state  of  statistical  control  as  quickly  as  possible. 

One  way  a  process  can  be  monitored  for  changes  is  through  the  use  of  a  control 
chart.  A  basic  Shewhart  control  chart  can  be  seen  in  Figure  1.  The  center  line  represents 
the  expected  in-control  value  of  the  quality  characteristic  being  measured.  Examples  of 
quality  characteristics  include  numerical  measurements  like  temperature  or  diameter,  or 
frequency  counts,  such  as  the  number  of  defects  inspected  per  item.  The  upper  control 
limit  (UCL)  and  lower  control  limit  (LCL)  are  also  shown.  The  control  limits  are  chosen 
in  such  a  way  that  if  the  process  is  statistically  in-control,  then  the  rate  at  which  the  chart 
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falsely  signals  is  controlled  at  some  desired  level.  In  Figure  1,  the  observations  are 
connected  by  a  line  to  show  how  the  observations  evolve  over  time. 


Figure  1:  Basic  Control  Chart 

Control  charts  are  closely  related  to  hypothesis  testing.  The  control  chart 
continually  tests  the  null  hypothesis  that  the  process  is  statistically  in-control.  If  the 
observations  stay  within  the  control  limits,  this  is  equivalent  to  failing  to  reject  the  null 
hypothesis.  Likewise,  a  point  plotted  outside  the  control  limits  is  equivalent  to  rejecting 
the  null  hypothesis.  For  example,  Figure  1  may  represent  some  process  where  the 
hypothesized  mean  of  the  quality  characteristic  being  measured  is  10  units.  The 
hypothesis  test  is 
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H0\jui=  10;  i  =  \,2,...,T 
10;  i  =  l,2,...,T. 

Looking  at  Figure  1,  one  can  fail  to  reject  the  null  hypothesis,  which  means  that 
there  is  not  enough  evidence  to  conclude  that  the  process  is  out-of-control.  Wackerly, 
Mendenhall,  and  Scheaffer  [22]  provide  an  excellent  overview  of  mathematical  statistics, 
including  hypothesis  testing. 

1.1.2  Shewhart  Control  Chart 

The  first  control  chart  used  in  this  thesis  is  the  Shewhart  control  chart,  named 
after  Walter  Shewhart,  who  developed  these  ideas  while  working  at  Bell  Labs  during  the 
1920’s.  A  general  model  used  to  construct  a  Shewhart  control  chart  is  as  follows.  Let  in 
denote  a  sample  statistic  measuring  some  quality  characteristic  of  interest.  The  mean  of 
m  is  /um  and  the  standard  deviation  of  m  is  am  .  The  center  line,  upper  control  limit,  and 
lower  control  limit  are 

UCL  =  pm+Lam 

CenterLine  =  jum  ( 1 . 1 ) 

LCL  =  -  L(J,n 

where  L  denotes  the  distance  of  the  control  limits  from  the  center  line  expressed  in 
standard  deviation  units.  A  control  chart  designed  according  to  these  principles  is  a 
Shewhart  control  chart,  and  using  L  =  3  creates  a  three-sigma  ( 3cr )  control  chart. 

Integrating  all  of  the  above  concepts  and  terminology,  an  SPC  process  monitored 
by  a  3cr  control  chart  can  be  described  as  follows.  In  phase  1,  process  data  are  gathered 
over  some  period  of  time  in  which  the  process  is  assumed  to  be  in-control.  The  data  are 
analyzed,  trial  control  limits  are  established  using  equation  (1.1),  and  process  engineers 
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then  use  the  control  chart  to  bring  the  process  into  a  state  of  statistical  control.  Once 
confident  the  process  is  in  control  and  the  actual  3 cr  control  limits  of  the  in-control 
process  have  been  calibrated,  then  phase  II  is  implemented.  During  phase  II,  the  process 
is  monitored  by  comparing  the  control  chart  statistic  computed  from  each  sample  drawn 
to  the  control  limits.  If  the  process  shifts  to  an  out-of-control  state,  the  control  chart  will 
likely  signal,  and  the  process  engineers  can  begin  their  search  for  any  assignable  causes 
that  may  have  caused  the  process  to  change.  Their  ultimate  goal  is  to  make  appropriate 
process  adjustments  to  get  the  process  back  into  a  state  of  statistical  control  as  quickly  as 
possible. 

One  way  to  evaluate  the  performance  of  a  control  chart  implemented  in  phase  II 
monitoring  is  to  look  at  the  Average  Run  Length  (ARL).  ARL  is  defined  as  the  number 
of  observations  plotted  on  the  control  chart  before  a  point  plots  outside  the  control  limits, 
signaling  an  out-of-control  condition.  For  a  Shewhart  control  chart  where  the 
observations  are  uncorrelated,  the  ARL  is  calculated  from 

ARL  =  —  (1.2) 

P 

where  p  is  the  probability  any  point  exceeds  the  control  limits.  For  a  3cr  control  chart 
where  the  mean  is  the  quality  characteristic  of  interest,  p  =  0.0027  is  the  probability  that 
a  point  falls  outside  the  control  limits.  The  average  run  length  of  the  chart  when  the 
process  is  in-control  (ARL0)  is 


ARL0 


- =  370. 

0.0027 
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The  interpretation  is  even  if  the  process  remains  in-control,  an  out-of-control  signal  will 
be  generated  on  average,  every  370  samples.  Obviously  a  large  ARL0  is  a  desirable 
property  of  a  control  chart. 

Shewhart  control  charts  are  designed  a  little  differently  depending  on  the  quality 
characteristic  being  measured.  The  two  main  types  of  Shewhart  control  charts  are  control 
charts  for  variables  and  attributes.  If  the  quality  characteristic  being  measured  can  be 
represented  on  a  numerical  scale,  such  as  temperature  or  length,  the  numerical  scale  is 

called  variable.  The  two  main  control  charts  for  variables  are  x  and  R  charts  (or  S 
charts)  used  to  monitor  the  mean  and  variance  of  numerical  quality  characteristics, 
respectively.  However,  not  all  quality  characteristics  are  represented  numerically. 
Perhaps  each  item  inspected  is  classified  as  confonning  (nondefective)  or  nonconforming 
(defective).  Quality  characteristics  of  this  type  are  called  attributes.  There  are  three 
main  Shewhart-type  control  charts  for  attributes: 

1.  /7-charts  for  binomial  counts.  Monitors  the  fraction  of  nonconforming  product 
produced  by  a  process. 

2.  c-charts  for  Poisson  counts.  Monitors  the  number  of  nonconformities  observed. 

3.  //-charts  for  Poisson  counts.  Monitors  the  average  number  of  nonconformities 
observed  per  unit. 

Again,  the  reader  is  referred  to  Montgomery  [9]  and  Ryan  [19]  for  more  detailed 
descriptions  of  SPC,  including  the  details  of  these  control  charts. 
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1.1.3  Cumulative  Sum  Control  Chart 


The  Shewhart  control  chart  is  easy  to  implement,  widely  used,  and  good  for 
detecting  large  changes  in  a  process.  It  is  particularly  useful  during  the  phase  I 
implementation  of  a  process,  where  the  process  is  not  finely  tuned  yet.  Shewhart  control 
charts  have  the  disadvantage  of  only  using  information  about  the  process  contained  in  the 
last  observation  while  ignoring  the  previous  sequence  of  observations.  Cumulative  sum 
(CUSUM)  control  charts  are  effective  alternatives  for  detecting  small  changes,  using 
information  from  the  current  and  previous  sequence  of  observations.  Therefore, 
CUSUMs  are  good  choices  to  implement  during  the  phase  II  monitoring  of  a  process 
when  the  in-control  process  parameters  are  assumed  to  be  known  or  sufficiently 
estimated.  Refer  to  Hawkins  and  dwell  [6]  for  a  complete  reference  on  CUSUM  charts, 
including  more  detail  regarding  the  theoretical  foundations  of  the  CUSUM  discussed 
below. 

The  tabular  CUSUM  plots  the  statistics 

Cf  =  max(0,  ±(y.-k)  +  Cf_l)  (1.3) 

where  Cf  is  the  cumulative  sum  at  time  i.  The  algorithm  cumulates  the  difference 

between  an  observed  value  yf  and  a  reference  value  k.  If  the  test  statistic  Cf  (Cf  ) 

exceeds  a  decision  internal  h+  (li  ) ,  the  chart  signals  that  an  increase  (decrease)  has 

been  detected,  suggesting  the  process  is  out-of-control.  The  values  of  /r  and  k 1  are 
selected  based  on  desired  ARL  properties. 

The  CUSUM  statistic  above  is  derived  from  the  sequential  probability  ratio  test 
(SPRT),  which  uses  observed  data  sequentially,  employing  each  observation  v,  as  it 
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becomes  available.  For  each  new  T,  the  SPRT  sequentially  tests  the  null  hypothesis 
H0 : 0(  =  0O  i  =  1,  2,  T  against  the  alterative  hypothesis  Ha :  0(  =  0a  i  =  1,  2, T  . 

0O  and  0<(  are  the  in-control  and  pre-specified  out-of-control  values  of  the  process 
parameters,  respectively.  Each  hypothesis  is  associated  with  a  PDF  /o(y,  |0o)  and 
fa(y,  I  0J  respectively. 

In  SPRT,  the  likelihood  ratio  is  given  by  Zf  =  ,  while  the  log- 

likelihood  ratio  is  defined  as  Cf  =  .  The  SPRT  operates  by  comparing 

Cf  to  a  decision  interval  A*  at  each  new  observation.  New  samples  are  collected  until 
Cf  >h±\  the  point  that  the  test  concludes  in  favor  of  Ha .  In  the  CUSUM  approach,  the 

null  hypothesis  of  in-control  is  never  accepted.  Hence,  the  test  is  restarted  each  time  the 
evidence  favors  the  null  hypothesis.  Thus,  samples  are  obtained  until  the  null  hypothesis 
is  rejected  in  favor  of  the  alternative  hypothesis. 

It  is  worth  noting  the  tabular  CUSUM  algorithm  uses  counters  (or  built-in  change- 
point  estimators)  that  suggest  when  a  process  shift  occurred.  For  example,  the  quantity 
N+  indicates  the  number  of  consecutive  periods  or  observations  that  the  CUSUM  Cf  has 
been  nonzero.  As  an  example,  suppose  that  the  control  chart  signaled  at  period  T  =  20, 
with  N+=  5.  This  would  indicate  that  five  consecutive  periods  have  elapsed  since  Cf 

rose  above  the  value  of  zero.  Thus,  it  is  likely  the  process  was  last  in-control  at  period  20 
-5  =  15,  and  that  the  shift  likely  occurred  between  periods  15  and  16.  After  reviewing 
the  basics  of  SPC  and  learning  how  a  3cr  control  chart  and  CUSUM  chart  are  used  to 
detect  changes  in  a  process,  the  problem  definition  can  now  be  defined. 


1.2  Problem  Definition 


It  was  discussed  above  that  a  control  chart  is  often  used  to  detect  a  change  in  a 
process.  Following  a  control  chart  signal,  knowledge  of  the  time  and  magnitude  of  the 
change  would  simplify  the  search  for  and  identification  of  the  assignable  cause.  In  this 
research,  emphasis  is  placed  on  count  processes  where  the  data  are  overdispersed,  which 
implies  that  the  observed  variance  is  greater  than  the  hypothesized  variance.  This  section 
will  first  discuss  overdispersion  in  general.  Then,  some  discussion  is  provided  that 
addresses  the  problem  of  estimating  a  process  change  point  where  the  observations  are 
overdispersed  count  data. 

1.2.1  Overdispersion 

Consider  a  modeling  scenario  where  the  quality  characteristic  of  interest  is  a 
count.  An  event  count  is  the  realization  of  a  non-negative  integer-valued  random  variable 
[4].  Examples  include  the  number  of  war  casualties  in  Iraq  per  day,  the  number  of 
defects  a  product  manufacturer  experiences  per  week,  or  the  number  of  accidents  a 
construction  company  observes  per  day.  It  is  well  known  that  the  Poisson  probability 
distribution  provides  a  reasonable  probability  distribution  to  model  count  data.  See 
Haight  [5]  for  a  comprehensive  reference  on  the  Poisson  distribution. 

Let  Yi  denote  a  random  variable  distributed  Poisson(  A,t ),  where  >  0  is  the  rate 

parameter  at  time  i  and  denotes  the  expected  value  of  T .  The  Poisson  probability  mass 
function  (PMF),  a  discrete  function  defined  for  y.  e  [0,  1,  2,  ...],  is 
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(1-4) 


W=y,)=f(y,  W=t-r~ 

u! 

where,  P{Yi  =  y.  \  A)  denotes  the  probability  that  the  random  variable  Y.  will  take  on  the 
value  y. ,  given  A. . 

Suppose  it  is  postulated  that  Ai  =  A0  for  all  i  =  1,2 .  Under  Poisson  model 
assumptions,  one  should  expect  that  E(Yi)  =  Var(Yi)  =  A0,  i.e.,  equidispersion,  for  all 
i  =  1,2 .  Unfortunately,  this  is  rarely  the  case  in  practice.  More  often  the  observed 
variance  of  a  random  sample  is  greater  than  the  theoretical  variance,  implying  that  the 
count  data  are  overdispersed.  If  Poisson  count  data  are  collected  over  time  in  discrete 
samples,  overdispersion  can  occur  as  a  result  of  uncontrollable  fluctuations  in  the  rate 
parameter  over  this  time  domain. 

For  example,  consider  Figure  2.  The  solid  line  represents  a  vector  of  Aj  's ,  where 

the  /'lh  element  is  the  value  of  the  rate  parameter  for  the  7th  sample.  This  vector  is 
generally  not  observable;  instead,  only  the  noisy  counts  can  be  observed  (represented  by 
circles  in  Figure  2). 
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Figure  2:  Overdispersed  Count  Data 

Under  the  postulated  model  that  Zt  =  Z0  for  all  i  =  ,  the  sample  mean  and 

sample  variance  of  the  100  observations  should  be  nearly  equal.  The  estimate  of  the 
mean  of  the  data  in  Figure  2  is  given  by  Zo  =  10.40  ,  whereas  the  variance  is  estimated  to 

be  =41.58.  This  overdispersion  is  a  direct  result  of  the  fluctuations  in  the  Zt '  s  over 

the  sampling  interval.  If  the  rate  parameter  remains  constant  over  the  sampling  interval, 
then  the  data  are  no  longer  overdispersed,  as  shown  in  Figure  3.  The  estimate  for  the 

mean  of  the  count  data  in  Figure  3  is  given  by  Zo  =  9.48 ,  whereas  the  variance  is 

estimated  to  be  0>  =  9.57 .  Clearly,  the  mean  of  the  data  is  much  closer  to  the  variance 
of  the  data  given  in  Figure  3. 
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Figure  3:  Equidispersed  Count  Data 

It  was  discussed  above  that  the  Poisson  probability  distribution  is  often  not  an 
adequate  means  to  model  overdispersed  count  data.  Again,  this  is  due  to  the  observed 
variance  of  a  random  sample  being  greater  than  the  theoretical  variance.  Thus,  one 
objective  of  this  research  is  to  derive  a  probability  distribution  that  will  adequately 
account  for  overdispersion  in  count  data  that  is  otherwise  distributed  as  Poisson.  In  this 
section  overdispersion  was  defined.  It  is  shown  in  Chapter  3  that  the  two-parameter 
gamma-Poisson  mixture  model  is  a  reasonable  model  for  overdispersed  count  data.  The 
next  section  will  discuss  the  problem  of  change-point  estimation  within  an  SPC  setting, 
where  the  count  data  are  assumed  to  be  overdispersed. 

1.2.2  Change  Point  Estimation 

Often  the  process  engineer  would  like  to  not  only  detect  a  change  in  the  process, 
but  also  estimate  when  the  change  actually  occurred.  Knowledge  of  the  time  and 
magnitude  of  the  change  would  simplify  the  search  for  and  identification  of  the 
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assignable  cause.  Having  an  estimate  for  the  change  point  can  reduce  the  costs 
associated  with  misdiagnosing  a  control  chart  signal.  These  costs  include  the  time  and 
resources  used  to  track  down  any  assignable  causes,  as  well  as  any  costs  associated  with 
unnecessary  process  adjustments.  Thus,  the  primary  focus  of  this  research  is  to  derive 
and  evaluate  maximum  likelihood  change-point  estimators  (MLEs)  for  the  time  of  a  step 
change  in  each  of  the  parameters  of  the  gamma-Poisson  mixture  model  following  a 
genuine  control  chart  signal.  After  a  control  chart  detects  the  process  change,  the 
engineer  can  very  quickly  isolate  when  the  true  change  may  have  occurred  by  using  the 
estimated  change  point  as  a  guide. 

1.3  Research  Objectives  and  Assumptions 

The  objectives  of  this  research  are  as  follows: 

1.  Derive  the  gamma-Poisson  mixture  model  that  will  be  used  throughout  this 
research  effort. 

2.  Derive  a  change-point  estimator  for  the  time  of  step  change  in  each  of  the 
parameters  of  the  gamma-Poisson  mixture  model  using  the  method  of 
maximum  likelihood  estimation. 

3.  Use  Monte  Carlo  simulation  to  evaluate  the  root  mean  square  error 
performance  of  the  change-point  estimators. 

4.  Present  and  evaluate  an  approach  based  on  the  likelihood  function  for 
estimating  a  confidence  set  for  the  process  change  point. 
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The  assumptions  made  in  this  research  are  as  follows: 


1.  This  research  effort  will  investigate  change-point  estimation  techniques  for 
the  two-parameter  gamma-Poisson  mixture  model.  In  the  first  case,  it  is 
assumed  that  the  process  experiences  a  step  change  in  the  location  parameter 
following  some  unknown  change  point  r  .  That  is,  the  process  starts  in  an  in¬ 
control  state  where  the  actual  in-control  mean  is  known  or  sufficiently 
estimated  from  a  phase  I  study.  After  some  unknown  point  in  time,  z ,  the 
process  mean  suddenly  shifts  to  some  unknown,  out-of-control  value  and 
remains  at  that  value  until  detected  and  corrected  by  the  process  engineer. 

2.  For  the  second  case,  the  assumption  is  that  the  process  experiences  a  step 
change  in  the  overdispersion  parameter  following  some  unknown  change 
point  z  .  That  is,  the  process  starts  in  an  in-control  state  where  the  actual  in¬ 
control  overdispersion  parameter  is  known  or  sufficiently  estimated  from  a 
phase  I  study.  After  some  unknown  point  in  time,  z ,  the  overdispersion 
parameter  suddenly  shifts  to  some  unknown,  out-of-control  value  and  remains 
at  that  value  until  detected  and  corrected  by  the  process  engineer. 

1.4  Thesis  Organization 

This  thesis  is  divided  into  five  chapters.  Chapter  1  presented  relevant  background 
information  on  SPC,  the  motivation  for  this  research,  and  the  problem  to  be  studied. 
Then,  the  research  objectives  and  assumptions  of  this  thesis  were  presented.  Chapter  2 
reviews  relevant  literature  concerning  change-point  estimation  within  an  SPC  context. 
Chapter  3  develops  the  mathematical  foundation  for  the  proposed  change-point 
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estimators.  Chapter  4  discusses  the  Monte  Carlo  simulation  study  that  was  used  to 
evaluate  the  performance  of  the  proposed  methodologies,  and  presents  the  results.  Last, 
Chapter  5  presents  the  conclusions  of  this  research  and  proposes  recommendations  for 
future  research. 
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2.  Literature  Review 


Pignatiello  [16]  points  out  that  very  little  attention  has  been  given  to  change-point 
estimation  in  the  quality  engineering  literature.  This  chapter  presents  relevant  change- 
point  estimation  research,  as  applied  within  an  SPC  context,  which  has  been  given 
attention  in  the  literature.  This  literature  shows  processes  that  experience  step  changes, 
linear  trends,  and  mono  tonic  changes  have  been  studied  for  normal,  Poisson,  and/or 
binomial  processes.  First,  this  chapter  starts  with  a  review  of  a  change-point  model 
designed  for  step  changes. 

2.1  Overview  of  a  Change-Point  Model 

The  change-point  model  used  in  this  research  is  attributed  to  Hinkley  [7]. 
Flinkley  proposed  a  change-point  model  for  step  changes  in  nonnal  process  means. 
Montgomery  [9]  states  that  the  change-point  procedure  is  very  good  at  estimating  the 
point  in  time  when  a  process  has  changed,  and  that  it  should  be  given  wider  attention  in 
process  monitoring  problems.  Montgomery  gives  an  introduction  to  the  change-point 
model,  under  nonnal  theory  assumptions,  applied  within  an  SPC  context.  Now  a  change- 
point  model  for  step  changes  is  defined. 

Consider  a  process  where  a  sustained  shift  has  occurred  in  a  parameter,  0  .  The 

change-point  model  is  defined  as 

T,  ~/o(©o)>  i  =  l,2,... ,t 
yi~fa(®a),i  =  t  +  U  +  2,...T 

where  the  in-control  distribution  is  /O(0O)  up  to  and  including  sample  r  ,  the  unknown 
change  point.  The  sustained  parameter  shift  occurs  between  r  and  z  + 1 ,  and 
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observations  are  sampled  from  the  out-of-control  distribution  fa  ( 0a )  from  time  r  +  1  to 

T.  Note,  0  can  represent  a  vector  of  parameters,  as  well  as  a  single  parameter.  Now  that 
the  basic  change-point  model  has  been  introduced,  some  relevant  literature  on  change- 
point  estimation,  within  an  SPC  context,  will  be  discussed. 

2.2  Change-Point  Estimation 

This  section  will  review  relevant  literature  on  change-point  estimation  considered 
in  the  quality  engineering  literature.  The  reader  is  referred  to  Basseville  and  Nikiforov 
[1]  for  a  more  general  coverage  of  change-point  detection  and  estimation  algorithms. 

Three  probability  models  often  used  in  SPC  applications  are  the  nonnal,  Poisson, 
and  binomial  distributions.  The  normal  distribution  is  probably  the  most  important 
distribution  in  the  theory  and  application  of  statistics.  For  instance,  normal  theory  is  the 
foundation  for  Shewhart  x  charts.  The  Poisson  probability  model  is  an  important 
discrete  distribution  used  in  SPC.  In  a  manufacturing  context,  for  example,  the  Poisson 
model  is  often  used  to  model  the  number  of  defects  (or  nonconformities)  that  occur 
within  an  inspection  unit.  Likewise,  the  binomial  probability  model  is  another  important 
discrete  distribution  used  in  SPC.  This  discrete  distribution  is  often  an  appropriate  model 
for  sampling  from  a  large  population,  where  parameter  p  represents  the  fraction  of 
defective  (or  nonconforming)  items  in  the  population.  The  relevant  change-point 
estimation  literature  concerning  these  distributions  is  discussed  next. 


17 


2.2.1  Normal  Case 


Samuel,  Pignatiello,  and  Calvin  [21]  consider  monitoring  a  process  with  a 
Shewhart  x  control  chart.  They  suggest  using  the  method  of  maximum  likelihood 
estimation  to  estimate  the  time  of  a  step  change  in  the  mean  of  the  process.  These 
authors  acknowledge  that  a  Shewhart  control  chart  can  signal  some  time  later  from  when 
a  change  actually  occurred  in  the  process  mean.  Thus,  they  argue  the  importance  of 
having  an  estimate  of  the  process  change  point.  They  derive  their  estimator  and  then 
evaluate  its  performance  using  Monte  Carlo  simulation.  They  show  that  applying  the 
MLE  following  a  genuine  Shewhart  control  chart  signal  provides  an  accurate  estimate  for 
the  change  point. 

Khoo  [8]  used  the  MLE  proposed  by  Samuel,  Pignatiello,  and  Calvin  [21]  to 
identify  the  time  of  a  permanent  shift  in  the  mean  of  a  process  being  monitored  by 
CUSUM  control  charts.  He  concluded  that  using  this  change-point  estimator  with  the 

CUSUM,  rather  than  the  Shewhart  x  control  chart,  offers  some  advantages.  Namely,  for 
small  shifts,  the  expected  run  length  is  lower  than  that  of  the  x  chart,  therefore  allowing 
process  improvement  to  be  carried  out  earlier. 

Pignatiello  and  Samuel  [16]  compare  the  MLE  of  the  time  of  a  step  change  in  a 
nonnal  process  mean  to  the  built-in  change-point  estimators  offered  by  the  CUSUM  and 
exponentially  weighted  moving  average  (EWM A)  control  charts.  Nishina  [10]  points 
out  that  the  CUSUM,  EWMA,  and  moving  average  charts  are  similar  with  respect  to 
estimating  a  process  change  point,  but  that  the  CUSUM  is  more  efficient  than  the  other 
two  in  change-point  estimation.  Using  Monte  Carlo  simulation,  Pignatiello  and  Samuel 
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[16]  show  that  the  performance  of  the  MLE  is  better  than  the  built-in  estimators  of  both 
the  EWMA  and  CUSUM  when  the  entire  range  of  step  change  magnitudes  is  considered. 
Not  only  do  they  consider  just  the  point  estimate,  but  they  also  propose  a  confidence  set 
for  the  time  of  the  process  change.  The  estimated  change  point  denotes  the  most  likely 
location  for  a  step  change  in  the  parameter  of  interest,  while  the  confidence  set  provides  a 
window  of  potential  locations  to  search  for  an  assignable  cause.  A  confidence  set  on  the 
process  change  point  is  an  advantage  of  using  the  MLE  as  opposed  to  the  CUSUM  or 
EWMA  change  point  estimators.  Note,  while  EWMA  control  charts  were  not  considered 
for  this  thesis,  an  interested  reader  may  refer  to  Roberts  [18]  or  Montgomery  [9]  for 
further  details  on  this  control  chart. 

Rather  than  assuming  the  out-of-control  process  behavior  is  adequately  modeled 
by  a  step  function,  Perry  and  Pignatiello  [12]  derive  the  MLE  for  the  time  of  linear  trend 
change  in  a  normal  process  mean.  They  compared  performance  results  between  their 
estimator  and  that  suggested  by  Samuel,  Pignatiello,  and  Calvin  [21]  following  genuine 

x  control  chart  signals  when  a  linear  trend  disturbance  was  present.  Perry  and 
Pignatiello  [12]  use  Monte  Carlo  simulation  to  evaluate  the  perfonnance  of  their 
estimator.  They  show  that  for  a  process  where  a  linear  trend  is  present,  the  MLE  of  the 
process  change  point  derived  for  linear  trends  outperforms  the  MLE  for  step  changes. 
They  also  evaluated  the  performance  of  the  proposed  estimator  when  the  distribution  of 
the  random  variable  being  monitored  is  not  exactly  normal.  It  was  shown  that  the 
estimator  perfonns  well  when  the  distribution  is  symmetric,  but  heavier-tailed  than  the 
nonnal.  Last,  they  presented  an  approach  using  the  likelihood  function  for  estimating  a 
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confidence  set.  This  would  provide  process  engineers  a  smaller  subset  of  observations  to 
investigate  for  an  assignable  cause  of  the  process  change. 

Park  and  Park  [11]  argue  that  the  x  chart  is  affected  by  changes  in  the  process 
variance  as  well  as  the  process  mean.  For  example,  many  of  the  authors  cited  above  will 
estimate  the  time  of  process  change  based  on  a  change  in  the  mean  alone.  Park  and  Park 
state  that  when  the  x  chart  issues  a  signal,  a  wrong  conclusion  could  be  reached  if  the 
change-point  estimator  considers  a  change  in  the  process  mean  only.  Thus,  they  apply  a 
joint  estimator  of  a  change  point  to  the  x  and  S  control  charts,  which  considers  the 
change  in  the  process  mean  and  variance  simultaneously.  Using  Monte  Carlo  simulation 
they  showed  that  the  proposed  estimator  perfonned  well  irrespective  of  subgroup  size 
over  various  changes  in  both  the  process  mean  and  variance. 

2.2.2  Poisson  Case 

Samuel  and  Pignatiello  [20]  derived  their  estimator  for  a  Poisson  process  under 
the  assumption  of  a  step  change  in  the  Poisson  rate  parameter  X .  They  applied  the 
proposed  estimator  after  a  signal  from  a  c- chart.  Using  Monte  Carlo  simulation  to 
evaluate  the  perfonnance  of  their  proposed  estimator,  they  concluded  overall  good 
accuracy  and  precision  perfonnance.  They  demonstrated  that  the  estimator  improves  as 
the  magnitude  of  the  step  change  increases. 

Perry,  Pignatiello,  and  Simpson  [15]  derived  and  evaluated  the  MLE  of  the 
process  change  point  using  the  change  likelihood  function  for  a  linear  trend  in  a  Poisson 
rate.  As  expected,  their  proposed  estimator  provides  good  overall  accuracy  and  precision 
performances  for  a  process  that  experiences  a  linear  trend  relative  to  the  MLE  derived  for 
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step  changes.  They  note  that  if  the  MLE  derived  for  step  changes  is  applied  to  processes 
that  experience  a  linear  trend,  the  estimator  will  tend  to  overestimate  the  true  change 
point. 

As  discussed  above,  many  different  change-point  methods  have  been  suggested. 
Most  of  these  methods  assume  the  type  of  change  in  a  process  parameter  is  known  before 
hand.  Perry,  Pignatiello,  and  Simpson  [14]  propose  a  maximum  likelihood  estimator  for 
the  time  of  change  of  a  Poisson  rate  parameter  without  requiring  exact  knowledge,  a 
priori,  regarding  the  exact  form  of  the  change  present.  They  evaluated  the  performance 
of  their  estimator  relative  to  that  achieved  by  the  estimator  suggested  by  Samuel  and 
Pignatiello  [20]  following  genuine  Poisson  CUSUM  control  chart  signals  and  across 
several  types  of  monotonic  behavior  in  the  rate  parameter.  These  monotonic  functions 
include  sudden  step,  linear  trend,  exponential,  and  log  functions.  Perry,  Pignatiello,  and 
Simpson  [14]  concluded  it  is  better  to  use  their  proposed  estimator  when  the  form  of  the 
change  present  is  only  known  to  be  monotonic.  Their  proposed  method  is  important 
since  process  engineers  rarely  know  the  type  of  change  before  hand. 

2.2.3  Binomial  Case 

The  discussions  above  on  change-point  estimation  have  focused  on  normal  and 
Poisson  processes  that  experience  step  changes,  linear  trends,  and  even  the  more  general 
monotonic  changes.  Regarding  change-point  estimation  for  binomial  processes, 
Pignatiello  and  Samuel  [17]  consider  a  step  change  in  p,  the  process  fraction 
nonconforming.  They  propose  an  estimator  using  the  method  of  maximum  likelihood 
estimation  to  estimate  when  the  change  in  p  began.  This  estimator  can  be  applied  after  a 
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p- chart  or  an  np-chart  has  detected  a  change.  Monte  Carlo  simulation  was  used  to  show 
that  the  estimator  provides  accurate  and  precise  estimates  of  the  time  of  a  step  change 
when  the  magnitude  of  the  change  in  p  is  at  least  10%  and  larger.  They  also  stress  the 
ease  of  using  this  method  in  a  spreadsheet  application. 

Perry  and  Pignatiello  [13]  also  considered  a  process  that  experiences  a  step 
change  in  the  process  fraction  nonconforming.  They  compared  the  performance  of  the 
MLE  proposed  by  Pignatiello  and  Samuel  [17]  to  built-in  change-point  estimators  from 
binomial  CUSUM  and  EWMA  control  charts.  They  concluded  that  the  MLE  of  the 
change  point  will  provide  an  overall  better  estimate  of  the  process  change  point  when 
compared  to  the  built-in  estimates  of  the  CUSUM  and  EWMA.  They  do  point  out  that 
the  only  exception  to  this  is  when  the  process  changes  to  the  pre-specified  magnitude  of 
change  used  to  design  the  CUSUM  or  EWMA  control  charts.  However,  since  these 
magnitudes  are  rarely  known  before  hand,  the  MLE  is  still  the  better  choice  overall. 

2.3  Conclusion 

This  chapter  started  with  a  review  of  a  change-point  model  for  step  changes. 
Next  this  chapter  reviewed  some  relevant  literature  on  change-point  estimation  as  applied 
within  an  SPC  context.  As  discussed,  processes  that  experience  step  changes,  linear 
trends,  and  monotonic  changes  in  general  have  been  studied  for  normal,  Poisson,  and/or 
binomial  processes.  However,  the  case  of  estimating  a  change  point  in  a  process 
parameter,  where  the  observations  are  overdispersed  counts,  has  currently  not  received 
much,  if  any,  attention  in  the  quality  engineering  literature.  Hence,  that  is  the  focus  of 
this  research.  The  methodology  of  this  research  will  be  developed  next. 


22 


3.  Methodology 


This  chapter  focuses  on  the  mathematical  models  used  for  this  research.  First,  the 
gamma-Poisson  mixture  model  is  derived.  It  is  the  probability  model  that  will  be  the 
basis  for  modeling  overdispersed  count  data  in  this  thesis.  Then,  two  different  step- 
change  scenarios  will  be  addressed.  The  change-point  model  for  a  process  that 
experiences  a  step  change  in  mean  is  defined,  followed  by  the  derivation  of  the  MLE  of 
the  process  change  point.  Next,  the  change-point  model  for  a  process  that  experiences  a 
step  change  in  the  overdispersion  parameter  of  the  underlying  probability  distribution  is 
defined.  Then,  the  corresponding  MLE  of  the  process  change  point  for  this  parameter 
will  be  derived.  Last,  the  method  of  building  a  confidence  set  to  capture  the  true  change 
point  is  discussed. 

3,1  Derivation  of  the  gamma-Poisson  Mixture  Model 

The  Poisson  distribution  was  defined  previously  in  equation  (1.4)  to  introduce  the 
idea  of  overdispersion.  For  completeness  of  deriving  the  gamma-Poisson  mixture  model, 
it  will  be  restated  here. 

Let  Yt  denote  a  random  variable  distributed  as  Poisson(  Xj ),  where  Xt  >  0  is  the 
rate  parameter  at  time  i  and  denotes  the  expected  value  of  Yj .  The  Poisson  probability 
mass  function  (PMF)  is  a  discrete  function  defined  for  yt  e  [0,  1,  2,  ...]  given  by 

P(Y,=y,)=ny,\K>=— '^r  O-i) 

yt ! 
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where,  P(Y.  =  y.  \  At)  denotes  the  probability  that  the  random  variable  Yj  will  take  on  the 
value  y, ,  given  2, . 

Suppose  that  it  is  postulated  that  At  =  X0  for  all  i  =  1,2,. ..,7  .  Then  under  Poisson 
model  assumptions,  one  should  expect  that  E(Yt)  =  Var(  Y  )  -  A0  ,  i.e.,  equidispersion,  for 
all  i  =  .  Unfortunately,  this  is  rarely  the  case  in  practice.  Instead  it  is  often  true 

that  the  observed  variance  of  a  random  sample  is  greater  than  the  theoretical  variance, 
implying  that  the  count  data  are  overdispersed.  If  Poisson  count  data  are  collected  over 
time  in  discrete  samples,  overdispersion  can  occur  as  a  result  of  uncontrollable 
fluctuations  in  the  rate  parameter  over  this  time  domain. 

It  is  reasonable  to  assume  the  variability  in  Aj  is  well  modeled  using  the  flexible, 

two-parameter  gamma  probability  distribution.  Let  2,  >  0  be  distributed  gamma  (a,/?) , 
where  /?  =  )),  oc  >  0  is  the  shape  parameter,  and  ft  >  0  is  the  scale  parameter.  The 
gamma  PDF  is  a  continuous  distribution  with  non-negative  support  defined  as 


f(Ai\a,p)  = 


A“-Xe  p 

r {a)p°  ‘ 


(3.2) 


If  Ai  denotes  a  random  variable  possessing  a  gamma  distribution  with  parameters 
a  and  /?  for  i  =  1,  2,  ...,  T  ,  then  it  is  easily  shown  that  E(A.)  =  v  and  Var(/ i(.)  =  -£ .  In 
equation  (3.2),  T(z)  denotes  the  gamma  function  where  for  a  non-negative  integer  z, 
T(z  + 1)  =  z ! ,  where  !  denotes  the  factorial  function. 
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The  joint  density  function,  using  the  parameterization  /?  =  £ ,  is  given  as 


f(yM  =  f(y,\WW 


r(y,  +  i)r(«) 

Thus,  the  marginal  distribution  of  y.  for  i  =  1,  2,  ...,  T  is 


f(y,)  =  \f(ylAl)d^ 

0 

-aAj 

r  dX 

'{  r(y,+l)T(a) 

r(i+f)(-^g)(frr(y,.+a) 
r(y,  +  i)r(«) 


n >/+a) 

f  v  ) 

(  a  ^ 

r(v;+l)T(a) 

^a  +  v  ) 

(3.3) 


(3.4) 


which  is  the  well-known  gamma-Poisson  mixture  parameterization  of  the  negative 
binomial  distribution  with  parameters  r  >  0 ,  where  r  =  a  and  p  e  [0, 1] ,  where  p  =  . 


This  model  is  employed  as  the  underlying  model  for  the  data  in  this  research.  The 
negative  binomial  model  in  equation  (3.4)  can  also  be  written  in  the  more  familiar  form 


as 


Ad  pY  p'  (3.5) 

where  the  constants,  =  !  ,  ,  ( 1  ,  and  p  =  —  .  For  the  remainder  of  this 

thesis,  the  expanded  fonn  shown  in  equation  (3.4)  will  be  used.  Note,  by  rewriting 
equation  (3.4)  as 
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/(»)  = 


r  (y,+a) 

f  1  > 

-a 

(  y.  \ 

V' 

T(a)(a  +  v)y‘ 

[l+aJ 

vh!/ 

and  calculating  lim  /'( v, )  yields 

<*-»  oo 


lim/(y;)  = 

ar-»oo 


e  vJ 


y,  • 


which  is  the  Poisson  probability  mass  function.  This  shows  that  the  Poisson  distribution 
is  a  limiting  form  of  the  negative  binomial  distribution  as  a  goes  to  infinity.  Therefore, 
a  can  be  interpreted  as  a  parameter  accounting  for  any  overdispersion  that  may  have 
occurred  in  the  data.  Hence,  a  is  referred  to  as  the  overdispersion  parameter.  Next,  the 
mean  and  variance  of  the  negative  binomial  distribution  are  calculated  from  the  moment 
generating  function. 

The  moment  generating  function  of  /(y,  )  is 


«,(<)  =  £(e»)  =  Ze!’/n,> 

y 


yt 


T  (Vi+a) 

f  v  7 

Yi 

(  a  ^ 

r(y,  +i)r(a) 

'ya  +  v  j 

(tt  +  Uy 

r  a  7 

a 

f  t  \( 

a+v-ev 

Ka  +v  j 

v  a  +  v  y 

a 


,  ,  •  (3-6) 

ya  +  v-ev  J 

The  first  moment  of  Y.  about  the  origin  (i.e.,  the  expected  value  or  mean)  is  calculated  as 


u[  =E(Yi)  =  mm(  0) 


=  v. 


(3.7) 
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Likewise,  the  second  moment  about  the  origin  of  Yj  is 


u2=E{Y2)  =  m(2\  0) 

=  v2+i  +  v.  (3.8) 


Thus,  the  variance  of  T  can  be  calculated  as 

Var(Yt)  =  E(Y?)-E(Yi)2 
=  (v2+^  +  v)-v2 

=  i+y.  (3.9) 

Note  that  if  lim(-^- +  v)  is  calculated,  Var(Yt)  =  v  =  E(Yj) .  Thus,  as  a  goes  to 

a->co 

infinity,  the  variance  approaches  the  theoretical  variance  of  a  Poisson  random  variable. 
Now  that  the  probability  model  for  this  research  has  been  derived,  the  change-point 
models  can  be  defined  and  the  proposed  estimators  derived. 


3.2  Change-Point  Model  for  a  Step  Change  in  Mean 

The  change-point  model  for  a  step  change  in  mean  is  parameterized  as  follows: 

»  ~  NB(a,  ^-),  i  =  l,2,...,r 

(3.10) 

yt~  NB(a,-^r),  i  =  r  +  l,r  +  2 


This  model  assumes  an  in-control  process  with  distribution  NB(ar,^-)  up  to  and 

including  a  point  z ,  the  unknown  change  point.  Assume  a ,  the  overdispersion 
parameter,  and  v0,  the  in-control  mean,  are  both  known  or  sufficiently  estimated. 

Between  times  z  and  z  + 1  the  process  shifts  to  a  new  level,  va  ^  v(l ,  and  the  process 
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continues  with  a  sustained  shift  in  mean  with  the  distribution  NB(a,-^-).  For  this 
research,  va  is  the  unknown  mean  after  the  process  experiences  the  step  change,  and  T 

denotes  the  time  of  the  control  chart  signal.  This  change-point  model  can  be  used  to 
derive  a  maximum  likelihood  estimator  for  the  process  change  point,  which  is  shown 
next. 


3.3  Derivation  of  the  MLE  When  va  is  Unknown 


The  proposed  MLE  of  the  change  point  is  denoted  as  zv .  Assuming  a  process 
change  point  at  z  ,  the  likelihood  function,  derived  from  equation  (3.4)  is 

m-y.  mo 

/= 1  i=r+ 1 

where  Ki  =  ny%rl)  • 

Assuming  there  has  been  a  single  change  point,  the  goal  is  to  estimate  z ,  the 
unknown  change  point.  The  MLE  of  z  is  the  value  of  z  that  maximizes  the  likelihood 
function  in  equation  (3.1 1),  or  equivalently  its  logarithm.  Taking  the  natural  logarithm  of 
equation  (3.11)  and  simplifying  (less  the  constant Kt)  yields  the  log-likelihood  function 


k  a+v0  ' 


va  |  y,)  =  t  yT  ln(^)  +  a ln(^)  +(T-z)  yTr  ln(^-)  +  a  ln(^-)  .  (3.12) 


Again,  the  objective  is  to  estimate  the  unknown  change  point  z  ;  however,  there 
are  two  unknowns  in  the  log-likelihood  function:  r  and  va .  An  estimate  for  va  is 
needed  to  maximize  the  log-likelihood  function.  If  the  value  of  z  was  known,  then  the 
value  of  va  that  maximizes  the  log-likelihood  function  is  given  by  the  solution  to 
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ai  yr.,(r-o[^-(^f](g+^)  (T_r)a  p 
dva  va  a  +  va 

Solving  this  in  terms  of  va  yields 

Va(f)  =  yT,  r 

where  yT  T  =  (T  -r) _y(. .  Making  this  substitution  into  the  log-likelihood  function 
of  equation  (3.12)  gives 

l(T\vi)  =  T  [yr  ln(^)  +  a  ln(^)]  +  (T  -  r)  [yr>r  ln(^-)  +  a  ln(^-)] .  (3.13) 

Now  the  log-likelihood  function  has  only  one  unknown  parameter,  the  change 
point  t.  Evaluating  equation  (3.13)  over  all  possible  integer  change  points 
t  e  [0,  —  1]  in  search  of  the  maximum  yields 

zv  =  argmax{/(z|y,)}  (3.14) 

0  <t<T 

where  zv  is  the  maximum  likelihood  estimate  of  the  last  observation  taken  from  the  in¬ 
control  process.  The  MLE  of  r  will  be  applied  to  data  obtained  following  a  control  chart 
signal.  To  make  this  concept  more  concrete,  consider  the  simulated  example  shown  next. 
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3. 3. 1  MLE  Example 


Consider  a  fabric  manufacturing  process.  The  quality  characteristic  being 
measured  is  the  number  of  defects  in  each  6’xlO’  section  of  fabric  produced.  The 
manufacturing  process  is  very  difficult  to  get  to  the  point  where  there  are  no  defects.  In 
fact  the  in-control  mean  number  of  defects  is  two  per  60  square  feet  of  fabric,  i.e.  v0  =  2. 
For  this  example,  the  overdispersion  parameter,  a  ,  is  assumed  to  be  10.  Calculating  the 
control  limits  for  the  3 a  control  chart  yields  UCL  =  6.64  and  LCL  =  0.  If  the  number  of 
defects  equals  7  or  greater,  then  the  fabric  will  be  considered  defective  and  will  be 
scrapped.  The  number  of  defects  per  unit  of  fabric  is  shown  in  Table  1. 

Table  1:  Observed  Count  at  Each  Observation  i 


Observation  Number  of  Observation  Number  of 
i  Defects  i  Defects 


1 

2 

15 

0 

2 

1 

16 

3 

3 

2 

17 

3 

4 

2 

18 

2 

5 

3 

19 

1 

6 

1 

20 

0 

7 

0 

21 

1 

8 

2 

22 

1 

9 

1 

23 

2 

10 

4 

24 

1 

11 

0 

25 

1 

12 

3 

26 

3 

13 

3 

27 

4 

14 

0 

28 

7 

Twenty  seven  units  of  fabric  were  inspected  and  determined  to  be  within  the 
control  limits  and  the  28th  unit  of  fabric  exceeded  the  UCL.  Thus,  the  3cr  control  chart 
signaled  at  time  T  =  28.  The  proposed  estimator  was  applied  after  the  control  chart 


30 


signaled.  In  applying  the  proposed  change-point  estimator,  each  potential  change  point,  t, 
needs  to  be  used  to  calculate  its  corresponding  log-likelihood  value  from  equation  (3.13). 

To  evaluate  equation  (3.13),  first  all  the  averages,  3v,  =  (r-<r'E'  ,+l  y, ,  and 

yt  =  (0_1Ei=:  y.  need  to  be  calculated.  Recall  in  this  example  that  T  =  28.  To  evaluate 
equation  (3.13)  for,  say,  potential  change  point  t  =  26,  calculate  yTt,  and  yt  as 

6  =  (2 r'YZx,y>  =  ¥  =  5.5  and  y2f>  =  (26)-'X“y,  =  =  1-6154. 

Bringing  these  values  and  the  other  known  parameters  into  equation  (3.13)  yields 
a  log-likelihood  value  of -142.82  for  t  =  26.  Evaluating  equation  (3.13)  over  all  t  will 

provide  the  value  of  t  that  achieves  the  maximum  log-likelihood  value,  which  is  rv  =  26  . 
All  log-likelihood  values  for  this  example  are  shown  in  Table  2. 

In  this  example  t=  26  gives  the  largest  log-likelihood  value.  That  is,  the  proposed 

estimator  calculates  the  change  point  at  zv  =  26  .  The  conclusion  is  that  observation  i  = 
26  is  the  last  observation  from  the  in-control  process  and  observation  i  =  27  is  the  first 
observation  from  the  out-of-control  process.  Process  engineers  would  be  instructed  to 
look  for  a  change  in  the  process  between  observations  26  and  27.  Note  that  the  simulated 
change  point  for  this  example  is  r  =  25.  In  this  example,  the  point  estimate  is  off  by 
one,  which  also  alludes  to  the  usefulness  of  confidence  sets  (covered  later)  to  capture  the 
true  change  point.  A  similar  example  was  done  by  Samuel,  Pignatiello,  and  Calvin  [21] 
for  a  step  change  in  the  normal  process  mean. 
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Table  2:  Log-Likelihood  Values  at  each  Potential  Change  Point 


Observation 

i 

Number 
of  Defects 

t 

Log  - 

Likelihood 

Value 

1 

2 

0 

-145.9449 

2 

1 

1 

-145.9423 

3 

2 

2 

-145.9807 

4 

2 

3 

-145.9794 

5 

3 

4 

-145.9780 

6 

1 

5 

-145.9296 

7 

0 

6 

-145.9747 

8 

2 

7 

-146.0133 

9 

1 

8 

-146.0133 

10 

4 

9 

-146.0024 

11 

0 

10 

-146.0016 

12 

3 

11 

-146.0012 

13 

3 

12 

-146.0133 

14 

0 

13 

-145.9992 

15 

0 

14 

-145.9986 

16 

3 

15 

-145.8752 

17 

3 

16 

-145.9460 

18 

2 

17 

-145.9947 

19 

1 

18 

-145.9928 

20 

0 

19 

-145.9245 

21 

1 

20 

-145.6326 

22 

1 

21 

-145.3570 

23 

2 

22 

-144.9582 

24 

1 

23 

-144.7837 

25 

1 

24 

-144.0714 

26 

3 

25 

-142.9806 

27 

4 

26 

-142.8196 

28 

7 

27 

-143.1652 

3.4  Change-Point  Model  for  a  Step  Change  in  the  Overdispersion  Parameter 

The  previous  section  defined  the  change-point  model  and  derived  the 
corresponding  estimator  for  the  case  where  a  step  change  occurred  in  the  mean  of  the 
process.  This  section  considers  the  case  where  the  process  experiences  a  step  change  in 
the  overdispersion  parameter,  a  .  The  change-point  model  is  parameterized  as  follows: 
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(3.15) 


v,.  ~  NB(a0,-^),  i  =  1,2,...,  z 
yt  ~  NB(aa,^),  i  =  t  + 1,  z  +  2, 


T. 


This  model  assumes  an  in-control  process  with  distribution  NB(a0,^;)  up  to 
and  including  the  unknown  change  point,  z  .  Assume  a0 ,  the  in-control  overdispersion 
parameter,  and  v ,  the  process  mean,  are  both  known.  Between  times  r  and  z  + 1  the 
process  experiences  a  shift  in  the  overdispersion  parameter  to  a  new  level,  aa  ^  aQ ,  and 

the  process  continues  with  a  sustained  shift  with  the  distribution  NB(afl,-^r).  For  this 
research,  aa  is  the  unknown  overdispersion  parameter  after  the  process  experiences  the 

step  change,  and  T  denotes  the  time  of  the  control  chart  signal.  This  change-point  model 
is  used  to  derive  a  maximum  likelihood  estimator  for  the  process  change  point. 


3.5  Derivation  of  the  MLE  When  aa  is  Unknown 

Assuming  a  process  change  point  at  z ,  the  likelihood  function  derived  from 

equation  (3.4)  is 


L(z,aa 


»>=n 


r(>’,+«0) 

r()i+i)ri«0) 


/  v  yv,  (  ao  \«o  TT- 

''O'0+l",  t«0+V't  X 


r (>’,+«n)  /  y  v v,  /  aa  'a, 

r(_yJ.+l)r(aa)  ^ aa+v '  k aa+v ) 


(3.16) 


Z=1  i=T+ 1 

The  goal  is  to  estimate  z ,  under  the  assumption  that  there  is  a  single,  unknown 
change  point.  The  MLE  of  z  is  the  value  of  z  that  maximizes  the  likelihood  function  in 
equation  (3.16),  or  equivalently  its  logarithm.  Taking  the  natural  logarithm  of  equation 
(3.16)  and  simplifying  (less  a  constant,  ln^y)  yields  the  log-likelihood  function 
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(3.17) 


l(j, aa  |  yt )  =  X [" ln ~ ln(«o  +  v)(y,  +  «o )  +  «o  In (a0 ) 

i= 1  L 

X  -  ln(a« + + «a) + «a  in(«j 


+ 


i=r+l 


Again,  the  objective  is  to  estimate  the  unknown  change  point,  z  ;  however,  there 
are  two  unknowns  in  the  log-likelihood  function  of  equation  (3.17):  z  and  aa .  A  value 
for  aa  is  needed  to  maximize  this  log-likelihood  function.  If  the  value  of  z  were  known, 
then  the  value  of  aa  that  maximizes  the  log-likelihood  function  is  given  by  the  solution 
to 


dl 

daa 


X  |V(yi  +  cca)-y/(cca) 

i  =  t+ 1 


( aa+v ) 


+  + 1 


=  0 


where  y/(x)  =  d'n^x)  • 


(3.18) 


Since  there  is  no  closed-form  solution  for  aa ,  Newton’s  method  will  be  used  to 
evaluate  equation  (3.18)  over  each  potential  change-point  value.  This  provides  an 
estimate  of  aa  for  each  potential  z.  Substituting  an  estimate,  oca(z) ,  for  each  an  in 
equation  (3.17)  yields 


Kt  I  yt )  =  X  ln  ~  ln(a°  +  v^y‘  +  a0 )  +  ao  ln(«o ) 


i= 1 
T 


+ 


Z  ln  ( E7K)1 )  “ ln  A  ^  +  yKy* +  ^  1 0)  ln(  t  0)) 


i=r+ 1 L 


(3.19) 


Now  the  log-likelihood  function  in  equation  (3.19)  has  only  one  unknown 
parameter,  the  change  point,  z.  Evaluating  equation  (3.19)  over  all  possible  integer 
change  points  t  e  [0, 1, ...,  T  - 1]  in  search  of  the  maximum  yields 


za  =  argmax  \l(t  |  y,)} 

0<t<T 


(3.20) 
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where  za  is  the  maximum  likelihood  estimate  of  the  last  observation  taken  from  the  in¬ 
control  process.  The  next  section  gives  a  short  summary  of  Newton’s  method  and  how  it 
was  used  to  estimate  aa  in  this  thesis. 


3.6  Newton’s  Method  for  Finding  the  MLE  of  aa 


Newton’s  method  is  a  derivative-based  algorithm  that  uses  a  linear  approximation 
to  find  roots  (zeros)  for  real-valued  functions.  Consider  finding  the  root  of  some 
differentiable  function  / .  Given  some  current  approximation,  xk ,  a  better 

approximation,  xk+l ,  may  be  found  using 


/(**) 


/'(**)' 


(3.21) 


Given  an  initial  guess,  x0,  chosen  near  the  root  of  the  function  /,  and  an 

appropriate  stopping  scheme,  equation  (3.2 1)  should  converge  at  the  root  of f  The  reader 
is  referred  to  Burden  and  Faires  [3]  for  more  detail  regarding  Newton’s  method.  Now, 
some  discussion  is  provided  on  the  application  of  Newton’s  method  to  the  problem 
presented  in  this  thesis. 

If  r  was  known,  Newton’s  method  could  be  used  to  solve  for  aa  in  equation 


(3.18).  The  estimate,  aa ,  at  the  k+ 1 st  iteration  can  be  calculated  explicitly  from 


'a,k+ 1 


'  a,k 


i  =  t+ 1 


T  r- 


%a,kyi+V 


i  =  t+\ 


0 \a,k&a,k+V )2 


(3.22) 
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where  aa  0  is  the  initial  guess,  <//(x)  =  digamma  function  =  Phi^x) ,  and  i//(l,x)  is  the 
trigamma  function.  For  this  thesis,  the  initial  guess  was  chosen  to  be  aa0  =  0.01 . 

Estimates  for  aa  can  be  obtained  at  each  potential  change  point,  t  e  [0,1,. ..7-1] , 


using  the  iterative  procedure  defined  in  equation  (3.22)  given  the  vector  of  observations 


y .  The  stopping  scheme  should  terminate  the  iterative  procedure  if  EL,*  -  EL,*+i  is 
sufficiently  small.  For  example,  if  the  square  root  of  the  sum  of  squares  of  the  difference 
of  the  last  two  iterations  (i.e.,  the  2-norm)  is  less  than,  say,  1.0  x  10"6.  The  estimator  from 
equation  (3.22)  is  then  used  for  the  evaluation  of  the  argument  of  equation  (3.20).  An 
example  of  Newton’s  method,  coded  in  MATLAB,  is  provided  in  Figure  4. 


a  k  =  0.01;  %Initial  guess 
tol  =  10A  (-6) ;  %Tolerance 
N  =  80;  %Number  of  iterations 

for  k  =  1 : N 

%Numerator  of  equation  (3.22) 

P=0 ; 

P  =  sum ( (psi ( 0 , y ( t+1 : T) +a_k)  -  psi(0,a_k)  -  ... 

( (y (t+1 : T) +a_k) . / (a_k  +  v) )  +  log (a_k/ (a_k+v) )  +  1) ) ; 

%Denominator  of  equation  (3.22) 

Q=0; 

Q  =  sum((psi(l,  y  (t+1 : T) +a_k)  -  psi(l,  a_k)  +  ... 

( (a_k.*y(t+l:T)  +  vA2 ) / ( ( (a_k+v) A2 ) *a_k) ) ) ) ; 

%Newton  steps 

a_kl  =  a_k  +  (-P/Q); 

error  =  norm(a  k  -  a  kl,  2); 

a  k  =  a  kl; 

if  error  <  tol 
break 

End 

if  a  k  <=  0  %Constraint  to  keep  a  k  positive 
a__k  =  10,000;  %Set  "large" 
break 

end 

end 


Figure  4:  Newton's  Method  Coded  in  MATLAB 
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3.6.1  Additional  Newton  Method  Constraint 

It  was  noticed  in  some  instances  of  Newton’s  method  that  aa  grows  very  large 

and  then  turns  negative  due  to  a  sign  change  in  the  numerator  of  equation  (3.22)  (denoted 
as  P  in  Figure  4).  However,  it  is  defined  that  a  >  0  .  Therefore,  the  additional  constraint, 
shown  in  the  code  of  Figure  4,  was  added  to  Newton’s  method  to  keep  the  estimates  of 
aa  positive.  Since  the  maximum  of  the  log-likelihood  function  of  (3.17)  in  these  cases  is 

achieved  for  very  large  aa ,  is  set  to  a  large  value  (i.e.,  =  10,000).  Figure  5  shows 

the  effect  of  increasing  a  over  the  log-likelihood  function,  as  well  as  the  numerator  (P) 
and  denominator  (Q)  of  equation  (3.22)  for  a  single  observation  y.  Notice  that  the 
denominator  (Q)  of  (3.22)  is  always  negative,  which  is  a  necessary  condition  for  (3.17)  to 
have  a  maximum. 


Figure  5:  Plots  Showing  the  Effect  of  Increasing  a  Over  the  Log-likelihood  Function,  and 
P  and  Q  from  MATLAB  code  in  Figure  4  for  Single  Observation  y.  v=5  and  y=5. 
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3.7  Confidence  Sets 


Estimating  r  only  gives  a  point  estimate.  This  research  now  constructs 
confidence  sets  on  the  process  change  points,  which  can  provide  the  process  engineer  a 
search  window  of  potential  change  points  that  covers  the  true  process  change  point  with  a 
given  level  of  confidence.  The  method  of  constructing  confidence  regions  on  parameter 
estimates  using  the  likelihood  function  is  attributed  to  Box  and  Cox  [2],  The  confidence 
set  is  defined  as  follows: 

CS  =  {Z :  In L(t)  >  In  L(z)  -  D)  (3.23) 

where  In  L(i)  is  the  maximum  of  the  log-likelihood  function  evaluated  over  all  potential 

change  points  t.  If  at  a  given  t,  the  value  of  In  L(t)  is  greater  than  lnZ(r)-Z) ,  where  D 
is  some  specified  constant,  then  t  is  included  in  the  confidence  set.  Box  and  Cox  suggest 
using  a  reference  value  D  =  \%2al  to  obtain  a  100(1  -a)%  confidence  region.  However, 

surface  plots  are  provided  in  Chapter  4  to  aid  the  user  in  selecting  a  value  for  D  that  will 
more  accurately  achieve  the  user’s  desired  level  of  confidence.  This  confidence  set 
approach  is  the  same  approach  used  by  Pignatiello  and  Samuel  [16]. 

3. 7. 1  Confidence  Set  Example 

To  make  equation  (3.23)  more  concrete,  consider  the  following  example.  A 
negative  binomial  process  has  in-control  parameters  NB(cr,-^r),  where  <2  =  10,  and 

v0  =  10 .  Following  the  50th  observation,  the  process  experiences  a  step  change  in  mean 

where  va  =14.  In  this  example,  zv ,  as  estimated  from  equation  (3.14),  yields  zv  =  47, 
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whereas  the  true  change  point  is  actually,  r  =  50.  Figure  6  is  a  plot  of  In  L(t)  versus  t 


using  a  3cr  control  chart  that  signaled  at  T  =  55.  Applying  the  confidence  set  estimator 
in  equation  (3.23)  with  D  =  1.5  yields  a  cardinality  of  6.  The  figure  shows  that  /  =  45 
through  t  =  50  are  in  the  confidence  set.  Note  the  true  change  point,  r  =  50,  is  included 
in  this  set.  A  confidence  set  will  give  process  engineers  a  window  of  potential  change 
points  that  they  can  use  to  search  for  an  assignable  cause.  Monte  Carlo  simulation  is 
used  to  evaluate  the  performance  of  these  confidence  sets  in  Chapter  4. 


Potential  Change  Point  (t) 


Figure  6:  Log-Likelihood  Plot 
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3.8  Conclusion 


This  chapter  started  by  deriving  the  gamma-Poisson  mixture  model.  This  is  the 
probability  model  that  was  the  basis  for  modeling  overdispersed  count  data  in  this 
research  effort.  In  addition,  two  different  step-change  scenarios  were  addressed.  First, 
the  change-point  model  for  a  process  that  experiences  a  step  change  in  mean  was  defined, 
followed  by  the  derivation  of  the  MLE  of  the  process  change  point.  Second,  the  change- 
point  model  for  a  process  that  experiences  a  step  change  in  the  overdispersion  parameter 
was  defined,  followed  by  the  derivation  of  the  MLE  for  this  case.  In  order  to  derive  the 
MLE  for  the  latter  case,  Newton’s  method  was  used  to  estimate  aa .  This  approach  was 

discussed,  along  with  the  additional  constraint  required  to  keep  the  estimates  of  aa 

positive.  Last,  the  method  of  building  a  confidence  set  to  capture  the  true  change  point 
was  discussed.  Chapter  4  uses  Monte  Carlo  simulation  to  evaluate  both  change-point 
estimators  proposed  in  this  chapter.  Likewise,  Monte  Carlo  simulation  is  used  to 
evaluate  the  perfonnance  of  the  confidence  set  estimators.  Chapter  4  concludes  with  an 
application  of  these  methods  to  a  real-world  problem. 
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4.  Results  and  Analysis 


In  Chapter  3,  a  change-point  estimator  was  derived  for  the  time  of  step  change  in 
each  of  the  parameters  of  the  gamma-Poisson  mixture  model  using  the  method  of 
maximum  likelihood  estimation.  This  chapter  discusses  the  Monte  Carlo  simulation 
experiments  used  to  evaluate  the  root  mean  square  error  perfonnance  of  the  proposed 

estimator,  rv ,  after  genuine  3cr  control  chart  signals.  Since  there  is  no  standardized 
negative  binomial  distribution,  this  thesis  investigates  the  12  combinations  of  the 
following  arbitrarily  chosen  in-control  parameter  values:  v0=  1 ,  5,  and  20,  and  <2=1,5, 
10,  and  50. 

Second,  this  chapter  discusses  the  Monte  Carlo  simulation  experiments  used  to 

evaluate  the  root  mean  square  error  performance  of  the  proposed  estimator,  r« ,  after 
genuine  CUSUM  control  chart  signals.  Since  there  is  no  standardized  negative  binomial 
distribution,  this  thesis  investigates  the  6  combinations  of  the  following  arbitrarily  chosen 
in-control  parameter  values:  v  =  5  and  20,  and  a0  =  1,5,  and  10. 

Monte  Carlo  simulation  is  also  used  to  evaluate  the  performance  of  the  confidence 
set  estimators  for  each  case  considered.  Last,  this  chapter  applies  the  methodologies 
developed  in  Chapter  3  to  Iraq  war  casualty  data  to  demonstrate  the  effectiveness  of  these 
techniques  on  a  real-world  problem. 

4,1  Monte  Carlo  Simulation  When  va  is  Unknown 

Monte  Carlo  simulation  is  used  to  evaluate  the  root  mean  square  error 
performance  of  the  change-point  estimator  following  a  genuine  signal  from  a  3cr  control 
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chart.  The  3cr  control  chart  was  chosen  due  to  its  wide  use  and  ease  of  implementation. 
Each  count  observed  will  be  compared  to  the  upper  and  lower  control  limits  of  the 
control  chart.  The  control  limits  are  calculated  as 


UCL  -  v0  +  3^v0  +  '-y- 
LCL=va-  3^7? 


(4.1) 


and  if  LCL  <  0,  then  set  LCL  =  0. 

The  process  change  point  was  simulated  to  occur  at  r  =  50.  Independent 
observations  were  drawn  from  a  negative  binomial  process  where  yi  ~~  NB(  for 

i  =  1,2,..., 50.  Following  the  50th  sample  collected,  observations  were  drawn  from  a 
negative  binomial  process,  where  jf~NB(a,^-)  for  i  =  r  + 1,  ...,  T,  where  T  is  the 


time  of  control  chart  signal.  Following  the  signal,  xv  was  calculated  from  equation 
(3.14).  This  procedure  was  repeated  a  total  of  A  =  10,000  times  for  each  value  of  va 

investigated.  The  overdispersion  parameter,  a ,  and  the  in-control  mean,  v0 ,  are  both 
assumed  known.  Note  that  a  remains  unchanged  after  the  step  change  occurs.  The  out- 
of-control  mean,  va ,  and  the  change  point,  z  ,  are  both  assumed  unknown. 


4.1.1  False  Alarms 

As  described  above,  each  simulated  observation  will  be  compared  to  the  upper 
and  lower  control  limits  of  the  3cr  control  chart.  If  the  control  chart  signals  at  an 
observation  T,  where  T  >  r ,  then  the  control  chart  signals  that  the  process  is  out-of¬ 
control.  However,  if  the  control  chart  signals  at  T,  where  T  <  r ,  then  it  is  treated  as  a 
false  alarm  since  the  chart  signaled  before  the  simulated  change  point.  In  the  simulation, 
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if  a  false  alarm  is  encountered  at  time  T,  the  control  chart  is  restarted  at  time  T+l,  and  the 
simulated  change  point  will  occur  as  scheduled.  This  is  consistent  with  how  a  false  alarm 
would  be  handled  in  practice,  i.e.,  if  a  false  alarm  is  encountered,  the  process  engineer 
would  declare  it  was  a  false  alann,  and  the  process  monitoring  is  continued  at  the  next 
observation.  Pignatiello  and  Samuel  [16]  handle  false  alarms  this  way. 

4.2  Performance  of  the  MLE  for  Step  Change  in  Mean 

The  performance  of  the  MLE,  rv ,  for  various  combinations  of  a  and  v0  are 
shown  in  the  tables  throughout  this  section.  Appendix  A  contains  all  the  results  from  the 
12  combinations  of  v0  and  a  studied  in  this  thesis.  Before  getting  to  the  specific  results 

shown  in  the  tables,  column  headings  of  the  tables  will  be  defined.  For  example,  see 
Table  3  below.  Recall  that  average  run  length  (ARL)  is  defined  as  the  expected  number 
of  observations  that  are  observed  before  the  control  chart  issues  a  signal.  Root  mean 
square  (RMS)  error  is  defined  as  the  square  root  of  mean  square  error  (MSE).  MSE,  a 
function  of  both  variance  and  bias,  is  defined  as 

MSE  =  +  Var(r)  (4.2) 

where  B(z)  =  E(t)  -  r  is  the  bias  of  the  estimator. 

Table  3  also  shows  rv ,  which  is  the  average  of  the  change-point  estimates 
obtained  over  the  A=  1 0.000  independent  simulation  runs.  The  last  column  of  the  table  is 

the  standard  error  of  zv  . 

Three  general  results  can  be  observed  from  the  data  in  this  section.  The  most 
obvious  is  for  a  given  v0  and  a  .  It  can  be  seen  that  by  looking  down  the  columns  of  any 
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individual  table,  say  Table  3,  that  as  the  magnitude  of  the  step  change  increases,  A  RL 

decreases,  rR[S(zv)  decreases,  and  zv  approaches  the  true  value  of  r  .  This  should  be 
intuitive  in  that  a  larger  step  change  in  mean  should  be  estimated  more  accurately  than  a 
smaller  step  change.  That  is,  the  bias  goes  to  zero  as  the  magnitude  of  the  step  change 
increases. 


Table3:  Performance  of  MLE  for  <2=1,  vQ  =  5,  r=50,  and  N  =  10,000 


K 

% 

Increase 

A&L 

R.RiS(zy) 

^  >| 

Std.error(  zv ) 

6 

20 

28.90 

37.49 

73.50 

0.292 

7 

40 

18.65 

21.99 

62.62 

0.180 

8 

60 

13.41 

14.78 

57.41 

0.128 

9 

80 

10.12 

11.23 

54.64 

0.102 

10 

100 

8.16 

9.07 

53.35 

0.084 

11 

120 

6.80 

7.96 

52.21 

0.076 

12 

140 

5.74 

7.08 

51.56 

0.069 

13 

160 

5.06 

6.40 

51.30 

0.063 

14 

180 

4.58 

5.96 

51.04 

0.059 

15 

200 

4.05 

5.93 

50.70 

0.059 

20 

300 

2.95 

4.52 

50.23 

0.045 

25 

400 

2.37 

4.33 

50.01 

0.043 

30 

500 

2.05 

4.06 

49.88 

0.041 

35 

600 

1.84 

4.15 

49.78 

0.042 

40 

700 

1.73 

3.28 

49.86 

0.033 

45 

800 

1.63 

3.55 

49.81 

0.035 

50 

900 

1.56 

3.57 

49.82 

0.036 

55 

1000 

1.47 

3.27 

49.80 

0.033 

The  second  general  result  that  can  be  seen  is  by  comparing  different  a's  for  a 
given  v0 .  For  example,  take  Table  3,  Table  6,  Table  7,  and  Table  8  and  compare  the 
summary  statistics  for  va  =35.  Table  4  summarizes  these  results  as  a  increases  for 
fixed  values  of  v0  =  5  and  vu  =  35  .  As  the  variance  of  the  observations  decreases  due  to 
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an  increasing  a ,  ARL  and  RMS(r,  )  both  decrease,  and  rv  approaches  the  true  value 
of  t  more  quickly.  Figure  7  visually  shows  the  effect  of  an  increasing  a  on  RMS. 


Table  4:  Effect  on  MLE  for  Increasing  a  . 

v0  =  5  and  va  =35.  r=50,  and  AMO, 000 


a 

A&L 

RSfc(rv) 

i< 

1 

1.84 

4.15 

49.78 

5 

1.08 

1.93 

49.81 

10 

1.02 

1.09 

49.92 

50 

1.00 

0.08 

50.00 

4.5  i  i  i  i  i  i  i  i  i 

4  - 

3.5  - 
3  - 

2.5  - 
if) 

2 

C  2 

1.5  - 

1  - 

0.5  - 

o _ 1 _ I _ 1 _ 1 _ 1 _ 1 _ 1 _ 1 _ 1 _ 

0  10  20  30  40  50  60  70  80  90  100 

alpha 

Figure  7:  Effect  on  RMS  due  to  an  Increasing  a  . 

v0  =  5  and  ua  =35.  r=50,  and  A=10,000 

A  similar  but  more  pronounced  example  of  this  second  result  can  be  seen  next. 
When  a  =  50,  the  estimator  approaches  the  true  change  point  much  more  quickly  than 
when  <2  =  1.  For  example,  when  a  =  50  and  v0  experiences  only  a  40%  change, 


45 


tv  =51.92.  This  is  in  contrast  to  when  a-  1  and  v0  experiences  a  40%  change.  In  this 

case,  Tv  =62.62.  These  results,  extracted  from  Table  3  and  Table  8,  show  the  effect  of 
alpha  on  the  MLE  and  are  summarized  in  Table  5  below.  Again,  the  conclusion  is  that  as 
the  overdispersion  parameter,  a ,  increases,  the  variance  of  the  observations  decrease, 
resulting  in  a  more  accurate  estimate  of  the  process  change  point. 


Table  5:  Comparison  of  MLE  for  <2  =  1  and  <2=50 . 
v0  =  5,  r=50,  and  A=10,000 


a  =  1 

l<  U 

a  ■ 

=  50 

Tv 

E,  =6 

73.50 

K 

=  6 

63.39 

& 

ii 

62.62 

va 

=  7 

51.92 

Va  =  10 

53.35 

K 

=  10 

49.77 

Table  6:  Performance  of  MLE  for  <2=5,  v0  =  5,  r=50,  and  N  =  10,000 


K 

% 

Increase 

A&L 

rRiS(zv) 

Tv 

Std.error(  xv ) 

6 

20 

39.01 

36.96 

70.98 

0.304 

7 

40 

19.45 

14.16 

55.62 

0.130 

8 

60 

11.39 

8.58 

51.83 

0.084 

9 

80 

7.36 

6.22 

50.67 

0.062 

10 

100 

5.28 

5.59 

50.03 

0.056 

11 

120 

4.12 

4.75 

49.92 

0.048 

12 

140 

3.31 

4.40 

49.71 

0.044 

13 

160 

2.78 

4.40 

49.57 

0.044 

14 

180 

2.41 

4.00 

49.61 

0.040 

15 

200 

2.15 

4.05 

49.50 

0.040 

20 

300 

1.49 

3.62 

49.55 

0.036 

25 

400 

1.25 

2.57 

49.70 

0.026 

30 

500 

1.14 

2.46 

49.71 

0.024 

35 

600 

1.08 

1.93 

49.81 

0.019 

40 

700 

1.06 

1.62 

49.84 

0.016 

45 

800 

1.03 

1.39 

49.89 

0.014 

50 

900 

1.02 

0.96 

49.93 

0.010 

55 

1000 

1.02 

1.13 

49.92 

0.011 
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Table  7:  Performance  of  MLE  for  a= 10,  v0  =  5,  r=50,  and  N  =  10,000 


K 

% 

Increase 

A&L 

rRiS(zy) 

l<  Is* 

Std.error(  rv ) 

6 

20 

51.07 

35.04 

68.49 

0.298 

7 

40 

22.41 

11.77 

53.68 

0.112 

8 

60 

11.93 

7.02 

50.98 

0.070 

9 

80 

7.28 

5.22 

50.26 

0.052 

10 

100 

4.91 

4.30 

49.89 

0.043 

11 

120 

3.66 

4.02 

49.69 

0.040 

12 

140 

2.86 

3.80 

49.64 

0.038 

13 

160 

2.36 

3.78 

49.53 

0.038 

14 

180 

2.07 

3.35 

49.60 

0.033 

15 

200 

1.79 

3.37 

49.54 

0.033 

20 

300 

1.26 

2.71 

49.65 

0.027 

25 

400 

1.11 

2.08 

49.78 

0.021 

30 

500 

1.04 

1.30 

49.89 

0.013 

35 

600 

1.02 

1.09 

49.92 

0.011 

40 

700 

1.01 

0.79 

49.95 

0.008 

45 

800 

1.00 

0.51 

49.97 

0.005 

50 

900 

1.00 

0.47 

49.98 

0.005 

55 

1000 

1.00 

0.70 

49.98 

0.007 

Table  8:  Performance  of  MLE  for  <2=50,  v0 

=  5,  T= 

=50,  and  A  =  10, 

K 

% 

Increase 

A&L 

rRB(zv) 

l<  Is* 

Std.error(  rv ) 

6 

20 

75.77 

30.75 

63.39 

0.277 

7 

40 

27.62 

9.54 

51.92 

0.094 

8 

60 

12.84 

5.63 

50.42 

0.056 

9 

80 

7.12 

4.27 

49.92 

0.043 

10 

100 

4.45 

3.54 

49.77 

0.035 

11 

120 

3.14 

3.36 

49.65 

0.033 

12 

140 

2.38 

3.26 

49.60 

0.032 

13 

160 

1.91 

3.11 

49.59 

0.031 

14 

180 

1.63 

3.01 

49.59 

0.030 

15 

200 

1.41 

2.67 

49.64 

0.027 

20 

300 

1.08 

1.61 

49.82 

0.016 

25 

400 

1.01 

0.95 

49.93 

0.010 

30 

500 

1.00 

0.33 

49.98 

0.003 

35 

600 

1.00 

0.08 

50.00 

0.001 

40 

700 

1.00 

0.05 

50.00 

0.001 

45 

800 

1.00 

0.02 

50.00 

0.000 

50 

900 

1.00 

0.00 

50.00 

0.000 

55 

1000 

1.00 

0.01 

50.00 

0.000 
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The  simulation  results  also  show  that  the  MLE  is  more  accurate  for  larger  in¬ 


control  values  of  v0  for  a  given  a  and  va .  For  example,  look  at  Table  8  along  with 
Table  10  and  Table  11.  It  can  be  seen  for  a  fixed  a  ,  and  a  larger  v0 ,  the  change  point  is 
estimated  much  more  accurately  for  a  given  percent  change  in  mean.  For  example,  a 
40%  increase  in  v0  =  1  is  va  =1.4,  which  gives  an  MLE  of  rv  =  1 15.84.  Contrast  this  to  a 

40%  increase  in  v0  =20  which  is  va  =28,  which  gives  an  MLE  of  zv  =52.56.  Note 

that  these  results  hold  due  to  the  fact  that  the  out-of-control  mean  is  specified  as  a 
percentage  increase  of  the  in-control  mean.  The  conclusion  here  is  that  a  larger  shift  in 
mean  will  result  in  a  more  accurate  estimate  of  the  true  change  point.  In  some  instances, 
in  practice,  smaller  shifts  (e.g.  10%  or  20%)  may  not  be  as  important  to  detect  as  large 
shifts.  Table  9  compares  the  MLE  for  different  in-control  values  of  v0  for  the  given 
percent  changes  shown. 


Table  9:  Comparison  of  MLE  for  a  =  50  for  v0  =  1,  5,  and  10. 

r=50,  and./V=10,000 


% 

Increase 

^o=1 

v0  =5 

%  =20 

20 

115.84 

63.39 

52.56 

40 

67.90 

51.92 

49.94 

100 

51.84 

49.77 

49.60 

400 

49.60 

49.93 

50.00 
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Table  10:  Performance  of  MLE  for  a=50,  v0  =  1,  r=50,  and  TV  =  10,000 


K 

% 

Increase 

A&L 

rWiS(zv) 

l<  Is* 

Std.error(rv ) 

1.2 

20 

117.43 

103.96 

115.84 

0.805 

1.4 

40 

64.79 

35.47 

67.90 

0.306 

1.6 

60 

39.11 

18.26 

57.18 

0.168 

1.8 

80 

25.67 

11.56 

53.28 

0.111 

2 

100 

17.80 

8.47 

51.84 

0.083 

2.2 

120 

13.02 

6.53 

51.05 

0.064 

2.4 

140 

9.96 

5.52 

50.57 

0.055 

2.6 

160 

7.81 

5.04 

50.29 

0.050 

2.8 

180 

6.35 

4.68 

50.06 

0.047 

3 

200 

5.26 

4.71 

49.87 

0.047 

4 

300 

2.68 

3.53 

49.65 

0.035 

5 

400 

1.81 

3.14 

49.60 

0.031 

6 

500 

1.44 

2.93 

49.65 

0.029 

7 

600 

1.23 

2.61 

49.67 

0.026 

8 

700 

1.13 

1.87 

49.78 

0.019 

9 

800 

1.07 

1.60 

49.83 

0.016 

10 

900 

1.04 

1.53 

49.86 

0.015 

11 

1000 

1.02 

1.29 

49.89 

0.013 

Table  11:  Performance  of  MLE  for  a= 50,  v0 

=  20,  r- 

=50,  and  A  =  10. 

K 

% 

Increase 

A&L 

R&S(r„) 

1  <  Is* 

Std.error(  zv ) 

24 

20 

28.56 

10.87 

52.56 

0.106 

28 

40 

7.76 

4.54 

49.94 

0.045 

32 

60 

3.33 

3.50 

49.61 

0.035 

36 

80 

1.96 

3.13 

49.56 

0.031 

40 

100 

1.45 

2.97 

49.60 

0.029 

44 

120 

1.21 

2.57 

49.66 

0.026 

48 

140 

1.10 

2.37 

49.73 

0.024 

52 

160 

1.04 

1.63 

49.84 

0.016 

56 

180 

1.02 

1.10 

49.91 

0.011 

60 

200 

1.01 

0.76 

49.94 

0.008 

80 

300 

1.00 

0.12 

49.99 

0.001 

100 

400 

1.00 

0.05 

50.00 

0.001 

120 

500 

1.00 

0.00 

50.00 

0.000 

140 

600 

1.00 

0.00 

50.00 

0.000 

160 

700 

1.00 

0.00 

50.00 

0.000 

180 

800 

1.00 

0.00 

50.00 

0.000 

200 

900 

1.00 

0.00 

50.00 

0.000 

220 

1000 

1.00 

0.00 

50.00 

0.000 
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General  results  indicate  that  as  the  magnitude  of  the  step  change  increases,  or  as 
the  variance  of  the  observations  is  lowered,  the  true  change  point  can  be  estimated  more 
accurately  in  a  process  that  experiences  a  step  change  in  mean.  The  above  results  show 
that  the  proposed  MLE  perfonns  well.  There  are  some  exceptions  for  small  step  changes 
with  highly  variable  data,  but  overall,  the  true  change  point  is  still  estimated  very 
accurately. 

The  importance  of  having  an  accurate  method  for  estimating  the  true  change  in  a 
process  should  now  be  more  obvious.  Look  at  Table  8  where  v0  =  5,  ua  =7,  a= 50,  and 

A&L=27.62 .  The  estimate  of  the  change  point  would  be  badly  biased  if  an  engineer 
estimated  the  time  of  the  process  change  at  the  time  the  control  chart  signaled.  It  is 
possible  that  the  process  engineer  could  incorrectly  diagnose  the  root  cause  of  the  change 
or  perhaps  not  even  discover  the  root  cause  at  all.  Instead  of  searching  for  causes  around 
the  correct  time  of  change,  a  process  engineer  might  examine  log  books  and  records 
associated  with  a  time  frame  that  was  almost  28  observations  after  the  process  actually 
changed. 

4.3  Cardinality  and  Coverage  Performances  of  Confidence  Set  Estimators  for 
Step  Change  in  Mean 

Set  cardinality  and  coverage  measures  are  used  to  evaluate  the  confidence  set 
estimators  for  the  nine  values  of  D,  where  D=  1,  1.25,  1.50,...,  3.  Recall  from  equation 
(3.23)  that  D  is  a  specified  constant  used  to  give  a  desired  confidence  and  cardinality  for 
each  percent  increase  in  v  relative  to  v0 .  The  change  point  was  evaluated  at  20% 

increases  in  the  mean  (relative  to  v0)  up  to  a  200%  increase.  Following  this,  the  change 
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point  was  evaluated  at  100%  increases  in  the  mean  up  to  a  1000%  increase.  Recall  that 
the  process  change  point  was  simulated  to  occur  at  r  =  50.  After  the  3cr  control  chart 
signaled  a  genuine  step  change,  the  confidence  set  estimator  was  applied.  The  cardinality 
of  the  confidence  set  was  recorded  as  well  as  whether  the  true  change  point  was  included 
in  that  set.  This  procedure  was  repeated  for  a  total  of  A  =  10,000  independent 

simulation  runs  for  each  percent  increase  in  v  relative  to  v0 .  The  average  cardinality 

was  computed  over  the  10,000  runs.  The  proportion  of  the  10,000  runs  that  included  the 
true  change  point  in  the  estimated  set,  i.e.,  coverage,  was  also  computed. 

Surface  plots  were  created  for  the  12  combinations  of  a  and  v0 .  These  plots  will 

aid  the  user  in  selecting  a  value  of  D  that  meets  his  or  her  desired  level  of  confidence. 
This  is  the  same  approach  used  by  Perry  and  Pignatiello  [12].  While  all  12  plots  are 
shown  in  Appendix  B,  four  of  the  surface  plots  are  shown  below. 

The  figures  show  that  more  coverage  can  be  obtained  for  a  given  percent  increase 
in  v  relative  to  v0  by  selecting  a  larger  reference  value  D.  The  tradeoff  is  that  the 

confidence  set  will  also  have  a  larger  cardinality.  For  example  in  Figure  8,  a  40%  change 
in  v  for  D  =  1.5  gives  a  coverage  of  0.39  and  cardinality  of  16,  whereas  D  =  2  gives  a 
coverage  of  0.64  and  cardinality  of  35.  For  the  case  where  D  =  2,  the  interpretation  is 
that  the  user  can  be  64%  confident  that  the  true  change  point  will  be  contained  in  the  set 
of  35  points. 
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Figure  8:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0 ,  and 

Reference  Value  D.  a  =  1,  v0=l,  r=50,  and  N  =  10,000. 


Figure  9:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  frorm/Q  ,  and 

Reference  Value  D.  a- 5,  v0=l,  r=50,  and  N  =  10,000. 
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Figure  10:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Increase  from  vQ ,  and  Reference 

Value  D.  a  =  10,  v0=l,  r=50,  and  N  =  10,000. 


Figure  11:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Increase  from  v0 ,  and  Reference 

Value  D.  a  =  50,  v0=l,  r=50,  and N  =  10,000. 
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4.3.1  Choice  of  z  for  Performance  Evaluation  of  the  Confidence  Set  Estimator 

For  this  research,  the  process  change  point  was  simulated  to  occur  at  r  =50. 
Results  in  Table  12  show  that  the  choice  of  z  in  this  simulation  study  has  little  impact  on 
the  coverage  probabilities.  Cardinality,  however,  appears  to  increase  slightly  for 
increases  in  z  ,  as  shown  in  Table  13.  While  different  values  of  z  could  have  been  used 
for  this  simulation  study,  overall  the  choice  of  r  =  50  seems  reasonable.  These  finding 
are  consistent  with  results  published  in  Pignatiello  and  Samuel  [16]. 


Table  12:  Effect  of  the  Change  in  z  on  Coverage  Probability. 

a  =  5,  v0  =20,  and  N  =  1 0,000 


%  Increase 

z  =  1 

o 
i r) 

II 

II 

o 

o 

20 

0.33 

0.28 

0.27 

80 

0.69 

0.64 

0.63 

160 

0.83 

0.81 

0.81 

300 

0.93 

0.91 

0.91 

Table  13:  Effect  of  the  Change  in  z  on  Average  Size  of  Confidence  Sets. 

a  =  5,  v0  =20,  and  N  =  \ 0,000 


%  Increase 

z  =  1 

o 

fn 

II 

II 

o 

o 

20 

5.42 

7.73 

8.26 

80 

2.28 

3.17 

3.47 

160 

1.59 

2.19 

2.42 

300 

1.29 

1.67 

1.72 

4.4  Monte  Carlo  Simulation  When  aa  is  Unknown 

Monte  Carlo  simulation  is  used  to  evaluate  the  root  mean  square  error 
performance  of  the  change-point  estimator  following  a  signal  from  a  CUSUM  control 
chart.  A  CUSUM  was  chosen  for  this  case  to  simply  add  variety  to  this  research  effort. 
The  derivation  of  this  CUSUM  is  in  section  4.5. 
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The  process  change  point  was  simulated  to  occur  at  r  =  50.  Independent 


observations  were  drawn  from  a  negative  binomial  process  where  y.  ~  NB(  oc0,-£^),  for 
i  =  1,2,..., 50.  Following  the  50th  sample  collected,  observations  were  drawn  from  a 
negative  binomial  process,  where  y.  ~NB(afl,-^j),  for  i  =  r  +  1,  ...,  T ,  and  T  denotes 

the  time  of  control  chart  signal.  Following  the  detection  signal,  ra  is  calculated  from 
equation  (3.20).  This  procedure  was  repeated  a  total  of  N  =  10,000  times  for  each  value 
of  aa  investigated.  Since  there  is  no  standardized  negative  binomial  distribution,  this 
thesis  investigated  the  6  combinations  of  the  following  arbitrarily  chosen  in-control 
parameter  values:  v  =  5  and  20,  and  a0=  1,  5,  and  10.  The  in-control  overdispersion 

parameter,  a0 ,  and  the  process  mean,  v ,  are  both  assumed  known,  where  v  remains 
constant.  The  out-of-control  overdispersion  parameter,  aa ,  and  the  change  point,  z  ,  are 
unknown.  False  alarms  were  handled  exactly  as  described  previously. 

4.5  CUSUM  Control  Chart  for  Detecting  a  Change  in  Overdispersion  Parameter 

For  this  research,  the  CUSUM  test  statistic  for  detecting  a  decrease  in  the 
overdispersion  parameter  is  given  by  C\  =  max{Cr,0} ,  where  C(  is  the  cumulative  sum 
at  time  i.  Evidence  of  a  decrease  in  the  overdispersion  parameter  is  indicated  by  a  signal 
occurring  if  C]  >  h  .  A  one-sided  tabular  CUSUM  was  designed  to  quickly  detect  a 
25%  decrease  in  a  relative  to  a0  with  the  in-control  average  run  length  (ARL0) 
calibrated  to  approximately  370. 
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The  CUSUM  control  chart  used  to  detect  a  decrease  from  a0  was  derived  from 

the  sequential  probability  ratio  test  (SPRT).  See  Hawkins  and  Olwell  [6]  for  details  of 
using  the  SPRT  to  derive  a  CUSUM  statistic. 

In  the  SPRT,  the  null  hypothesis  H0 :  a=a0,  i  =  1,  2, T  is  tested  against 


Ha:  a=au,  i  =  1,  2,  T  .  In  SPRT  ,  the  likelihood  ratio  is  given  by 


Z:  =■ 
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where  the  value  a0  is  the  in-control  value  of  a  ,  the  value  a*  is  the  pre-specified  out-of¬ 
control  value  of  a  that  one  wishes  to  detect,  and  v  is  the  process  mean,  assumed 
constant. 

To  derive  the  CUSUM  statistic,  take  the  natural  log  of  equation  (4.3)  and  simplify 
to  obtain 
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Rewriting  this  expression  gives 
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where  the  value  of  C0  is  taken  to  be  zero.  Equation  (4.5)  is  set  to  zero  whenever  C,  <  0 . 
The  SPRT  operates  by  comparing  Cr  to  the  decision  interval  h~  at  each  new  observation, 
if  c:  >  h~ ,  then  the  test  concludes  in  favor  of  Ha . 
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Note  that  the  familiar  tabular  form  of  the  CUSUM  statistic  is  not  ready  available 
since  the  observed  count  is  contained  in  the  gamma  function.  However,  since  this  is  a 


SPRT,  it  still  remains  a  CUSUM. 


4.6  Performance  of  the  MLE  for  Step  Change  in  the  Overdispersion  Parameter 

Three  general  results  can  be  observed  from  the  data  in  this  section.  The  most 
obvious  is  for  a  given  a0  and  v  .  It  can  be  seen  that  by  looking  down  the  columns  of  any 

of  these  tables,  say  Table  14,  that  as  the  magnitude  of  step  change  increases,  ARL  and 

RSlS(ra)  both  decrease,  and  za  approaches  the  true  value  of  r  .  The  conclusion  is  that 
a  larger  magnitude  of  step  change  in  a  will  result  in  a  more  accurate  estimate  of  the 
change  point,  i.e.,  the  bias  goes  to  zero  as  the  magnitude  of  the  step  change  increases. 


Table  14:  Performance  of  MLE  for  v  =  5,  «0=10,  z  =  50,  and  iV=10,000 


%  Decrease 

aBl  Rktef.) 

Za 

Std.error(  za ) 

9 

10 

246.97 

327.69 

276.17 

2.371 

8 

20 

175.21 

226.82 

205.20 

1.654 

7 

30 

124.48 

153.54 

155.15 

1.119 

6 

40 

85.02 

98.24 

115.80 

0.730 

5 

50 

59.27 

63.04 

90.23 

0.485 

4 

60 

40.55 

39.48 

72.81 

0.322 

3 

70 

26.70 

23.51 

60.72 

0.209 

2 

80 

16.92 

13.65 

53.77 

0.131 

1 

90 

9.26 

7.75 

50.18 

0.078 

The  second  general  result  that  can  be  seen  is  by  comparing  different  v's  for  a 
given  a0  and  given  percent  decreases.  For  example,  take  Table  14  and  Table  19  and 

compare  the  performance  of  za  for  v  =  5  and  v  =  20 .  Table  15  shows  that  when 
v  =  20 ,  the  estimator  approaches  the  true  change  point  more  quickly  than  when  v  =  5  . 
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The  conclusion  is  that  a  larger  v  will  experience  a  larger  magnitude  of  change  compared 
to  a  smaller  v  for  a  given  a0  and  given  percent  decrease.  Again,  the  bias  goes  to  zero  as 

the  magnitude  of  the  step  change  increases.  Note  that  the  results  in  this  section  hold  due 
to  the  fact  that  the  out-of-control  overdispersion  parameter  is  specified  as  a  percentage 
decrease  of  the  in-control  parameter. 

Table  15:  Comparison  of  MLE  for  v  =  5  and  v=20  for  Various  Decreases  in  a  Relative  to 

a0.  a0  =  10,  r=50,  and  N  =  10,000 


20% 

50% 

80% 

V 

Decrease 

Decrease 

Decrease 

5 

205.20 

90.23 

53.77 

20 

144.19 

64.21 

50.93 

The  third  general  result  from  this  section  is  that  the  MLE  is  more  accurate  for 
smaller  in-control  values  of  a0  for  a  given  mean  and  percent  decrease  in  a  relative  to 

a0 .  For  example,  look  at  Table  17,  Table  18,  and  Table  19  below.  It  can  be  seen  for  a 
fixed  v  and  given  percent  decreases  in  a  (relative  to  a0 ),  the  change  point  is  estimated 
more  accurately  at  a  lower  a0 .  For  example,  consider  a  50%  decrease  from  a  =  10  to 

a  =  5 ,  which  gives  an  MLE  of  ra  =  64.21 .  Contrast  this  to  a  50%  decrease  from  a  =  1 , 
to  a  =  0.5 .  This  case  gives  an  MLE  of  ra  =  56.3 1 .  The  conclusion  is  that  a  smaller  a0 
causes  higher  variance  in  the  observations,  y,  for  a  given  mean  and  percent  decrease.  It 

appears  that  this  increased  magnitude  in  variance  allows  Ta  to  be  estimated  more 
accurately.  Note  that  the  variance  of y,  as  shown  in  equation  (3.9),  can  be  made  larger  in 
two  ways.  The  first  way  to  increase  the  variance  is  by  increasing  v  for  a  given  a  .  The 
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second  way  the  variance  can  be  made  larger  is  to  decrease  a  for  a  given  v ,  which  is 


what  was  just  discussed.  These  results  are  summarized  in  Table  16. 


Table  16:  Comparison  of  MLE  for  v  =  20  for  a0  =  1,  5,  and  10. 


r=50,  and  N  =  10,000 


%  Decrease 

a0=l 

II 

o 

e 

£ 

O 

II 

o 

20% 

109.05 

126.01 

144.19 

50% 

56.31 

60.60 

64.21 

80% 

49.51 

50.51 

50.93 

Table  17:  Performance  of  MLE  for  v  =  20,  a0  =  1,  r  =  50,  and  A=T0,000 


aa 

%  Decrease 

A&L 

R^LS(ra) 

53 

|<  Is* 

Std.error(  ra ) 

0.9 

10 

151.83 

198.58 

184.50 

1.461 

0.8 

20 

76.51 

89.16 

109.05 

0.668 

0.7 

30 

43.77 

46.07 

77.89 

0.367 

0.6 

40 

27.73 

25.38 

63.51 

0.215 

0.5 

50 

18.87 

15.70 

56.31 

0.144 

0.4 

60 

13.30 

10.27 

52.61 

0.099 

0.3 

70 

9.46 

7.22 

50.53 

0.072 

0.2 

80 

6.77 

5.68 

49.51 

0.057 

0.1 

90 

4.76 

4.57 

48.98 

0.045 

Table  18: 

Performance  of  MLE  for  v  =  20,  a0 

=  5,  r  =  50,  andA=10,000 

au 

%  Decrease 

A&L 

RfilS(rff) 

53 

|<  Is* 

Std.error(  r« ) 

4.5 

10 

177.80 

235.03 

209.76 

1.724 

4 

20 

93.53 

114.08 

126.01 

0.851 

3.5 

30 

55.23 

61.16 

88.53 

0.475 

3 

40 

35.67 

36.15 

70.57 

0.297 

2.5 

50 

23.46 

21.45 

60.60 

0.186 

2 

60 

15.53 

13.51 

54.73 

0.127 

1.5 

70 

10.52 

9.17 

52.00 

0.090 

1 

80 

6.77 

6.23 

50.51 

0.062 

0.5 

90 

3.99 

4.11 

49.86 

0.041 
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Table  19:  Performance  of  MLE  for  v  =  20,  a0  =  10,  r  =  50,  and  AMO, 000 


au 

%  Decrease 

A&L 

R^S(ra) 

53 

|<  Is* 

Std.error(  r« ) 

9 

10 

192.08 

252.41 

223.78 

1.831 

8 

20 

112.28 

140.27 

144.19 

1.039 

7 

30 

66.85 

76.22 

99.57 

0.579 

6 

40 

43.35 

45.26 

77.37 

0.360 

5 

50 

28.19 

26.67 

64.21 

0.226 

4 

60 

19.11 

16.79 

57.04 

0.152 

3 

70 

12.62 

10.87 

52.88 

0.105 

2 

80 

8.02 

7.47 

50.93 

0.074 

1 

90 

4.44 

4.76 

50.01 

0.048 

General  results  indicate  that  as  the  magnitude  of  the  step  change  increases,  the 
true  change  point  can  be  estimated  more  accurately.  On  average,  the  results  show  that 
the  proposed  MLE  performs  well.  There  are  some  exceptions  as  were  shown,  but  overall, 
the  true  change  point  is  still  estimated  very  accurately.  See  Appendix  C  for  all  cases 
investigated. 

4.7  Cardinality  and  Coverage  Performances  of  Confidence  Set  Estimators  for 
Step  Change  in  Overdispersion  Parameter 

Set  cardinality  and  coverage  measures  are  used  to  evaluate  the  confidence  set 
estimators  for  the  nine  values  of  D,  where  D=\,  1.25,  1.50,...,  3.  Recall  from  equation 
(3.23)  that  D  is  a  specified  constant  used  to  provide  a  desired  confidence  level  and 
cardinality  for  each  percent  decrease  in  a  relative  to  a0 .  The  change  point  was 

evaluated  at  10%  decreases  from  a0  up  to  a  90%  decrease.  As  before,  the  process 

change  point  was  simulated  to  occur  at  r  =50.  Following  a  genuine  CUSUM  control 
chart  signal,  the  confidence  set  estimator  was  applied.  The  cardinality  of  the  confidence 
set  was  recorded  as  well  as  whether  the  true  change  point  was  included  in  that  set.  This 
procedure  was  repeated  for  a  total  of  N  =  10,000  independent  simulation  runs  for  each 
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percent  increase  in  a  (relative  to  a0 ).  The  average  cardinality  was  computed  over  the 

10,000  runs.  The  proportion  of  the  10,000  runs  that  included  the  true  change  point  in  the 
estimated  set,  i.e.,  coverage,  was  also  computed. 

Surface  plots  were  created  for  the  six  combinations  of  v  and  a0  looked  at  in  this 

thesis.  The  plots  of  various  D  values  will  aid  the  user  in  selecting  a  D  value  that  meets 
his  or  her  desired  level  of  confidence.  While  all  six  plots  are  shown  in  Appendix  D,  two 
of  the  surface  plots  are  shown  below. 

As  before,  the  figures  show  that  more  coverage  can  be  obtained  for  a  given 
percent  decrease  in  a  relative  to  a0  by  selecting  a  larger  reference  value  D.  The 
tradeoff  is  that  the  confidence  set  will  also  have  a  larger  cardinality.  For  example  in 
Figure  12,  a  50%  decrease  from  a0  for  D  =  1.5  gives  a  coverage  of  0.48  and  cardinality 

of  27,  whereas  D  =  2  gives  a  coverage  of  0.64,  but  a  cardinality  of  42.  For  the  case  where 
D  =  2,  the  interpretation  is  that  the  user  can  be  64%  confident  that  the  true  change  point 
will  be  contained  in  the  set  of  42  points. 
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Cardinality 


Figure  12:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=5,  a0= 5,  r=50,  and  TV  =  10,000. 


Figure  13:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=20,  a0=  5,  r=50,  and iV  =  10,000. 
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4. 7.1  Choice  of  z  for  Performance  Evaluation  of  the  Confidence  Set  Estimator 


For  this  research,  the  process  change  point  was  simulated  to  occur  at  t  =  50. 
Results  in  Table  20  show  that  the  choice  of  z  in  this  simulation  study  has  little  impact  on 
the  coverage  probabilities.  Cardinality,  however,  appears  to  increase  slightly  for 
increases  in  r  ,  as  shown  in  Table  21.  While  different  values  of  z  could  have  been  used 
for  this  simulation  study,  overall  the  choice  of  z  =  50  seems  reasonable.  These  finding 
are  consistent  with  the  results  published  in  Pignatiello  and  Samuel  [16]. 


Table  20:  Effect  of  the  Change  in  z  on  Coverage  Probability. 

a  =  5,  v0  =20,  and  N  =  1 0,000 


%  Decrease 

z  =  1 

O 

fT) 

II 

II 

o 

o 

10 

0.10 

0.09 

0.09 

50 

0.45 

0.38 

0.39 

90 

0.71 

0.65 

0.65 

Table  21:  Effect  of  the  Change  in  z  on  Average  Size  of  Confidence  Sets. 

a  =  5,  v0  =20,  and  N  =  1 0,000 


%  Decrease 

z  =  1 

O 

in 

II 

o 

o 

II 

10 

17.39 

18.62 

19.17 

50 

7.31 

10.57 

11.47 

90 

1.91 

2.83 

2.99 
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4.8  Change-Point  Analysis  Applied  to  Iraq  War  Coalition  Casualty  Data 

The  data  and  results  of  this  thesis  have  been  obtained  from  simulation 
experiments.  Hence,  it  begs  the  question,  “How  are  these  techniques  applied  to  a  real-life 
data  set?”  To  answer  this  question,  a  data  set  containing  Iraq  war  coalition  casualties  was 
obtained  from  the  Iraq  Coalition  Casualty  Count  website  at  http://icasualties.org/oii7.  A 
realistic  scenario  using  that  data  is  as  follows. 

The  US  commander  in  Iraq  has  tasked  the  strategic  assessments  team  to  monitor 
coalition  casualties  and  to  report  if  there  are  any  noticeable  changes  from  the  average 
daily  coalition  casualties.  The  assessment  team  decided  they  would  use  a  3 a  control 
chart  to  monitor  for  changes,  and  that  they  would  also  like  to  apply  change -point 
estimation  techniques  after  the  control  chart  signals  a  genuine  change.  Not  knowing 
exactly  what  kind  of  change  might  occur;  they  decided  to  employ  a  step  change  model. 

Because  daily  casualties  are  data  counts,  and  it  is  suspected  that  the  data  would  be 
overdispersed  due  to  the  variability  associated  with  war,  the  gamma-Poisson  mixture 
model  was  used  as  the  underlying  probability  distribution  of  the  counts.  The  assessment 
team  was  to  start  monitoring  on  07/19/2003.  To  get  a  baseline,  the  equivalent  of  the 
phase  I  SPC  calibration  process,  the  assessment  team  took  prior  data  from  a  100-day 
period  from  04/10/2003  -  07/18/2003.  Recall  that  the  gamma-Poisson  mixture  model,  as 
parameterized  in  this  thesis,  is  the  negative  binomial  probability  distribution  with 
parameters  r  =  a  and  p  =  .  After  some  careful  analysis,  the  team  estimated  the  in¬ 

control  parameters  as  &  =  2.21  and  Vo  =  1.25  .  They  established  upper  and  lower  control 
limits  as  UCL  =  5.4  and  LCL  =  0  respectively.  Any  daily  coalition  casualty  count  of  6  or 
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greater  would  cause  the  control  chart  to  signal.  They  were  now  ready  to  use  the  control 
chart  to  monitor  daily  coalition  casualties. 

Figure  14  shows  the  results  of  the  strategic  assessment  team’s  change-point 
analysis.  The  team  started  monitoring  daily  coalition  casualties  starting  on  07/19/2003 
(day  1).  The  control  chart  detected  a  change  on  10/26/2003,  or  day  100  from  when  the 
daily  monitoring  was  started.  After  the  control  chart  signaled,  they  applied  the  MLE 
from  equation  (3.14)  and  estimated  the  change  point  as  day  99,  as  shown  graphically  in 
Figure  15.  Recall  that  the  change  point  is  the  last  point  from  the  in-control  process.  In 
this  case,  the  control  chart  signal  and  the  estimated  change  point  are  in  exact  agreement. 

The  log-likelihood  plot  of  Figure  15  is  important  in  the  fact  that  it  can  be  used  to 
search  for  assignable  causes  on  other  likely  days.  Even  though  the  MLE  of  the  change 
point  was  day  99,  it  can  be  seen  from  Figure  15  that  potentially  any  of  the  change  points 
from,  say,  days  83-99  can  be  investigated  for  assignable  causes.  As  an  example  of  the 
usefulness  of  the  confidence  set  estimator,  the  surface  plot  of  Figure  8  is  used  as  a  guide 
to  narrow  down  the  pool  of  potential  change  points  even  further.  By  setting  D  =  1.75 
from  equation  (3.23),  the  assessment  team  can  narrow  their  search  window  for  assignable 
causes  to  about  12  days  and  still  achieve  almost  a  90%  confidence. 

After  the  control  chart  signal,  there  are  two  conclusions  that  can  be  made.  Either 
this  signal  is  a  false  alann  and  nothing  is  out  of  the  ordinary,  or  that  this  is  indeed  a 
genuine  signal  and  it  warrants  further  investigation.  The  team  decided  to  investigate 
further  and  they  discovered  that  the  Muslim  holy  month  of  Ramadan  begins  on 
10/27/2003,  or  day  101  from  when  daily  monitoring  began.  Their  conclusion  was  that 
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this  is  not  a  false  alarm  and  they  briefed  their  results  to  the  US  commander  in  charge  so 
that  appropriate  actions  could  be  taken. 

Figure  16  is  insightful  in  that  it  shows  what  happened  after  the  change  was 
detected  in  Figure  14.  In  Figure  14,  even  though  a  change  was  detected,  the  daily 
casualty  does  not  seem  to  be  too  far  out  of  the  ordinary.  But  as  Figure  16  shows,  the 
control  chart  indeed  detected  a  significant  change  that  is  clearly  shown  by  the  dramatic 
increase  in  casualties  during  the  month  of  Ramadan.  Flad  the  proposed  method  actually 
been  used,  tighter  security  measures  could  have  been  implemented  following  the  control 
chart  signal,  and  perhaps  coalition  lives  could  have  been  saved. 
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Figure  14:  Control  Chart,  Iraq  Coalition  Casualties  07/19/03  -  10/26/03 
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Figure  15:  Log-Likelihood  Plot  of  Iraq  War  Casualty  Data 
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In-control  values:  In-control  for  99  days 

alpha  =  2.21  07/1 9/03  -  1 0/25/03 

mean  =  1 .25  (days  1  -99) 


Step  change  in  mean 
10/26/03  (day  100) 

Ramadan  starts  10/27/03  (day  101) 
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Figure  16:  Coalition  Casualties  07/19/03  -  12/31/03 
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This  example  was  based  off  of  real-world  data.  A  3cr  control  chart  was  used  to 
monitor  a  process  where  the  data  was  overdispersed  count  data.  This  example  shows  the 
power  of  using  change-point  methods  in  conjunction  with  a  control  chart  to  monitor  a 
process.  The  control  chart  will  allow  a  change  to  be  detected,  while  change-point 
estimation  techniques  can  narrow  down  the  search  for  an  assignable  cause.  In  this 
example,  however,  the  estimated  change  point  and  control  chart  signal  coincided  largely 
due  to  a  spike  in  fatality  counts.  While  most  of  the  literature  deals  with  examples  from  a 
manufacturing  context,  this  example  shows  that  these  methods  can  be  applied  to  non- 
traditional,  non-manufacturing  processes,  including  military  processes. 

4.9  Conclusion 

This  chapter  employed  Monte  Carlo  simulation  to  evaluate  the  perfonnance  of 
each  of  the  MLEs  derived  in  Chapter  3.  The  first  case  considered  was  when  a  count 
process  experienced  a  step  change  in  the  mean.  Monte  Carlo  simulation  was  used  to 
evaluate  the  root  mean  square  error  performance  of  the  change-point  estimator  following 
a  signal  from  a  3cr  control  chart.  General  results  indicate  that  as  the  magnitude  of  the 
step  change  increases,  or  as  the  variance  of  the  observations  is  lowered,  the  true  change 
point  can  be  estimated  more  easily  in  a  process  that  experiences  a  step  change  in  mean. 

The  second  case  considered  was  when  a  count  process  experienced  a  step  change 
in  the  overdispersion  parameter.  Monte  Carlo  simulation  was  used  to  evaluate  the  root 
mean  square  error  performance  of  this  change-point  estimator  following  a  signal  from  a 
CUSUM  control  chart.  General  results  indicate  that  as  the  magnitude  of  the  step  change 
increases,  the  true  change  point  can  be  estimated  more  accurately. 
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The  performance  of  the  confidence  set  estimator  for  each  case  was  also 
investigated.  The  surface  plots  generated  will  aid  a  user  in  selecting  a  value  for  D  that 
meets  his  or  her  desired  level  of  confidence.  Last,  the  methodologies  developed  in 
Chapter  3  were  successfully  applied  to  Iraq  war  coalition  casualty  data  to  demonstrate  the 
effectiveness  of  these  techniques  on  a  real-world  problem.  Next,  Chapter  5  will 
summarize  this  research  effort  and  make  recommendations  for  future  research. 


69 


5.  Conclusions  and  Future  Research 


5.1  Summary  and  Conclusions 

If  a  process  engineer  could  determine  the  actual  time  of  a  process  change,  their 
time  and  resources  could  be  better  used  to  investigate  the  assignable  causes  that 
contributed  to  the  change.  A  quick  identification  of  the  assignable  causes  would  lead  to 
an  improvement  of  the  process  sooner.  Estimating  the  true  change  point  is  also  important 
due  to  the  costs  associated  with  misdiagnosing  a  control  chart  signal.  These  costs  include 
the  time  and  resources  used  to  track  down  any  assignable  causes,  as  well  as  any  costs 
associated  with  unnecessary  process  adjustments.  Therefore,  the  primary  focus  of  this 
research  was  to  derive  and  evaluate  maximum  likelihood  estimators  for  the  time  of  a  step 
change  in  each  of  the  parameters  of  the  gamma-Poisson  mixture  model,  following  a 
genuine  control  chart  signal. 

To  motivate  a  methodology  for  this  research,  a  review  of  the  relevant  literature 
was  presented.  First  was  a  review  of  a  step  change,  change-point  model.  Then,  the 
relevant  literature  on  change-point  estimation  was  discussed.  That  is,  processes  that 
experience  step  changes,  linear  trends,  and  monotonic  changes  have  been  studied  for 
normal,  Poisson,  and/or  binomial  processes.  Many  cases  of  change-point  estimation  are 
considered  in  the  literature;  however,  the  case  of  estimating  a  change  point,  where  the 
observations  are  overdispersed  count  data,  has  not  received  much,  if  any,  attention  in  the 
quality  engineering  literature. 

The  methodology  focused  on  the  mathematical  models  employed  in  this  research. 
The  gamma-Poisson  mixture  model  was  derived  and  employed  as  the  probability  model 
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for  modeling  overdispersed  count  data  in  this  thesis.  Then,  two  different  step-change 
scenarios  were  addressed.  First,  the  change-point  model  for  a  process  that  experiences  a 
step  change  in  mean  was  defined,  followed  by  the  derivation  of  the  MLE  of  the  process 
change  point.  Second,  the  change-point  model  for  a  process  that  experiences  a  step 
change  in  the  overdispersion  parameter  of  the  underlying  probability  distribution  was 
defined.  Then,  the  corresponding  MLE  of  the  process  change  point  was  derived  for  this 
case.  Last,  the  method  of  building  a  confidence  set  to  capture  the  true  process  change 
point  was  discussed.  Monte  Carlo  simulation  was  used  to  evaluate  each  of  these 
proposed  estimators.  Likewise,  Monte  Carlo  simulation  was  used  to  evaluate  the 
performance  the  proposed  confidence  set  estimators. 

The  first  case  considered  a  step  change  in  the  mean  of  a  process.  The  conclusions 
reached  were  as  the  magnitude  of  the  step  change  increased,  or  as  the  variance  of  the 
observations  was  lowered,  the  true  change  point  was  estimated  more  accurately.  Then  set 
cardinality  and  coverage  measures  were  used  to  evaluate  the  confidence  set  estimators. 
The  given  surface  plots  show  more  coverage  can  be  obtained  for  a  given  percent  increase 
in  v  from  v0  by  selecting  a  larger  reference  value  D,  but  the  tradeoff  is  that  the 

confidence  set  will  also  have  a  larger  cardinality. 

Second,  Monte  Carlo  simulation  was  used  to  evaluate  the  perfonnance  of  the 
proposed  estimator  after  a  CUSUM  control  chart  detected  a  step  change  in  the 
overdispersion  parameter.  Results  indicate  that  as  the  magnitude  of  the  step  change 
increased,  the  true  change  point  could  be  detected  more  accurately.  The  set  cardinality 
and  coverage  measures  were  also  used  to  evaluate  the  confidence  set  estimators  in  this 
case,  with  the  similar  findings  as  the  step  change  in  mean  case. 
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To  demonstrate  the  effectiveness  of  these  techniques  on  a  real-world  problem,  the 
methodology  of  Chapter  3  was  applied  to  Iraq  war  coalition  casualty  data.  A  3cr  control 
chart  was  used  to  detect  a  change  in  daily  casualties,  and  then  the  change-point  estimator 
was  applied.  It  was  shown  that  an  increase  in  daily  coalition  casualties  coincided  with 
the  start  of  the  Muslim  holy  month  of  Ramadan.  The  significance  of  the  example  is 
while  most  of  the  literature  deals  with  examples  from  manufacturing  settings,  this 
example  shows  these  methods  can  be  applied  to  a  non-traditional,  non-manufacturing 
process. 

5.2  Future  Research 

Although  a  general  methodology  has  been  developed  for  estimating  a  step  change 
in  both  the  mean  and  overdispersion  parameters  of  an  overdispersed  count  process,  there 
are  areas  of  future  research.  For  instance,  in  this  thesis,  the  two-parameter  gainma- 
Poisson  mixture  model  was  used  to  model  overdispersed  count  data.  It  was  assumed  that 
both  parameters  were  known  for  the  in-control  process,  but  that  the  process  went  out-of¬ 
control  due  to  a  step  change  in  only  one  of  the  parameters.  While,  these  assumptions 
allowed  a  foundation  to  be  laid,  other  cases  can  be  considered. 

The  case  where  neither  the  in-control  nor  out-of  control  parameter  values  are 
known  would  be  a  very  valuable  addition  to  this  research.  While  any  established  process 
would  have  known  or  sufficiently  estimated  in-control  parameter  values,  it  is  possible  a 
new  process  could  benefit  from  this  kind  of  analysis.  For  example,  at  the  start  of  the  Iraq 
war,  there  may  not  have  been  any  historical  data.  Certainly  estimating  a  change  in  a 
process  where  none  of  the  parameters  are  known,  a  priori,  is  important. 
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In  addition,  the  assumption  of  a  sustained  step  change  can  be  too  simplistic  in 
some  instances.  The  techniques  in  this  thesis  lend  themselves  well  to  a  process  where  it 
is  assumed  that  the  type  of  change  can  be  accurately  modeled  by  a  step  function. 
However,  if  the  type  of  change  is  not  known  before  hand,  then  the  procedure  proposed  by 
Perry,  Pignatiello,  and  Simpson  [14]  may  be  a  better  alternative. 

In  this  thesis,  it  was  assumed  that  the  variability  of  the  Poisson  rate  parameter,  A  , 
was  modeled  more  accurately  by  allowing  it  to  be  modeled  as  the  flexible,  two-parameter 
gamma  probability  distribution.  However,  it  may  be  the  case  that  a  set  of  covariates 
(predictors)  do  well  in  explaining  some  of  the  variability  in  A .  A  change-point  model 
that  incorporates  covariates  may  have  the  potential  to  model  a  wide  range  of  process 
behaviors.  For  example,  step  changes,  linear  trends,  and  cyclical  changes. 

Last,  this  research  proposes  that  the  negative  binomial  distribution  is  an  adequate 
model  to  use  for  all  count  data  processes  since  the  Poisson  distribution  is  just  a  limiting 
form  of  the  negative  binomial  distribution.  However,  the  one -parameter  Poisson 
distribution  may  be  preferred  in  cases  where  it  can  be  used.  Hence,  it  may  be  worth 
characterizing  for  exactly  what  values  of  various  combinations  of  the  negative  binomial 
parameters  that  the  Poisson  distribution  is  an  adequate  model  to  use  in  place  of  the 
negative  binomial. 
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Appendix  A:  Simulation  Data  for  Step  Change  in  Mean 


Table  22:  Performance  of  MLE  for  a  =  1,  v0=l,  z  =  50,  and  N  =  10,000 


va 

%  Increase 

A&L 

R  kfc(zv) 

|<  is* 

Std.error(  zv ) 

1.2 

20 

38.20 

50.42 

82.57 

0.385 

1.4 

40 

25.66 

30.89 

69.03 

0.243 

1.6 

60 

18.45 

20.97 

61.82 

0.173 

1.8 

80 

14.22 

15.46 

57.93 

0.133 

2 

100 

11.47 

12.00 

55.43 

0.107 

2.2 

120 

9.44 

10.10 

53.90 

0.093 

2.4 

140 

8.03 

8.70 

53.06 

0.082 

2.6 

160 

7.16 

7.76 

52.38 

0.074 

2.8 

180 

6.23 

7.08 

51.87 

0.068 

3 

200 

5.60 

6.61 

51.46 

0.065 

4 

300 

3.86 

5.27 

50.59 

0.052 

5 

400 

2.98 

4.72 

50.22 

0.047 

6 

500 

2.54 

4.10 

50.12 

0.041 

7 

600 

2.25 

3.99 

49.97 

0.040 

8 

700 

2.06 

3.60 

49.98 

0.036 

9 

800 

1.87 

3.83 

49.83 

0.038 

10 

900 

1.77 

3.49 

49.84 

0.035 

11 

1000 

1.70 

3.23 

49.86 

0.032 

Table  23: 

Performance  of  MLE  for  <2  =  1,  v0  =  5, 

r  =  50,  and  A  =  10,000 

K 

%  Increase 

A&L 

rRb(tv) 

|<  is* 

Std.error(  zv ) 

6 

20 

28.90 

37.49 

73.50 

0.292 

7 

40 

18.65 

21.99 

62.62 

0.180 

8 

60 

13.41 

14.78 

57.41 

0.128 

9 

80 

10.12 

11.23 

54.64 

0.102 

10 

100 

8.16 

9.07 

53.35 

0.084 

11 

120 

6.80 

7.96 

52.21 

0.076 

12 

140 

5.74 

7.08 

51.56 

0.069 

13 

160 

5.06 

6.40 

51.30 

0.063 

14 

180 

4.58 

5.96 

51.04 

0.059 

15 

200 

4.05 

5.93 

50.70 

0.059 

20 

300 

2.95 

4.52 

50.23 

0.045 

25 

400 

2.37 

4.33 

50.01 

0.043 

30 

500 

2.05 

4.06 

49.88 

0.041 

35 

600 

1.84 

4.15 

49.78 

0.042 

40 

700 

1.73 

3.28 

49.86 

0.033 

45 

800 

1.63 

3.55 

49.81 

0.035 

50 

900 

1.56 

3.57 

49.82 

0.036 

55 

1000 

1.47 

3.27 

49.80 

0.033 
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Table  24:  Performance  of  MLE  for  a  =  1,  vo=20,  z  =  50,  and  N  =  10,000 


K 

%  Increase 

A&L 

Rkfc(rv) 

|< 

Std.error(r,. ) 

24 

20 

28.52 

36.64 

72.91 

0.286 

28 

40 

17.86 

20.68 

61.66 

0.171 

32 

60 

12.52 

14.00 

56.84 

0.122 

36 

80 

9.48 

10.45 

54.21 

0.096 

40 

100 

7.57 

8.72 

52.80 

0.083 

44 

120 

6.19 

7.55 

51.93 

0.073 

48 

140 

5.52 

6.60 

51.57 

0.064 

52 

160 

4.85 

6.10 

51.19 

0.060 

56 

180 

4.29 

5.91 

50.85 

0.059 

60 

200 

3.90 

5.37 

50.66 

0.053 

80 

300 

2.76 

4.80 

50.09 

0.048 

100 

400 

2.26 

4.55 

49.89 

0.046 

120 

500 

1.99 

3.82 

49.93 

0.038 

140 

600 

1.79 

3.76 

49.83 

0.038 

160 

700 

1.67 

3.40 

49.87 

0.034 

180 

800 

1.58 

3.90 

49.72 

0.039 

200 

900 

1.50 

2.98 

49.84 

0.030 

220 

1000 

1.45 

2.90 

49.84 

0.029 

Table  25: 

Performance  of  MLE  for  a  =  5,  v0  =  1 , 

r  =  50,  and  A  =  10,000 

K 

%  Increase 

aSl 

R  fciS(Tv) 

|<  Is* 

Std.error(  r,- ) 

1.2 

20 

59.63 

69.91 

95.41 

0.532 

1.4 

40 

35.53 

31.62 

68.00 

0.260 

1.6 

60 

23.10 

17.99 

58.16 

0.160 

1.8 

80 

16.20 

11.95 

54.42 

0.111 

2 

100 

12.29 

9.13 

52.63 

0.088 

2.2 

120 

9.33 

7.42 

51.61 

0.072 

2.4 

140 

7.60 

6.37 

50.97 

0.063 

2.6 

160 

6.40 

5.68 

50.58 

0.057 

2.8 

180 

5.40 

5.25 

50.27 

0.053 

3 

200 

4.66 

5.09 

50.04 

0.051 

4 

300 

2.73 

4.10 

49.65 

0.041 

5 

400 

1.99 

3.79 

49.57 

0.038 

6 

500 

1.63 

3.47 

49.58 

0.035 

7 

600 

1.43 

3.23 

49.61 

0.032 

8 

700 

1.31 

2.76 

49.66 

0.027 

9 

800 

1.23 

2.70 

49.68 

0.027 

10 

900 

1.17 

2.32 

49.74 

0.023 

11 

1000 

1.12 

2.01 

49.78 

0.020 
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Table  26:  Performance  of  MLE  for  a  =  5,  v0  =  5,  r  =  50,  and  TV  =  10,000 


K 

%  Increase 

A&L 

R  fcte(Tv) 

'i  >1 

Std.error(r,- ) 

6 

20 

39.01 

36.96 

70.98 

0.304 

7 

40 

19.45 

14.16 

55.62 

0.130 

8 

60 

11.39 

8.58 

51.83 

0.084 

9 

80 

7.36 

6.22 

50.67 

0.062 

10 

100 

5.28 

5.59 

50.03 

0.056 

11 

120 

4.12 

4.75 

49.92 

0.048 

12 

140 

3.31 

4.40 

49.71 

0.044 

13 

160 

2.78 

4.40 

49.57 

0.044 

14 

180 

2.41 

4.00 

49.61 

0.040 

15 

200 

2.15 

4.05 

49.50 

0.040 

20 

300 

1.49 

3.62 

49.55 

0.036 

25 

400 

1.25 

2.57 

49.70 

0.026 

30 

500 

1.14 

2.46 

49.71 

0.024 

35 

600 

1.08 

1.93 

49.81 

0.019 

40 

700 

1.06 

1.62 

49.84 

0.016 

45 

800 

1.03 

1.39 

49.89 

0.014 

50 

900 

1.02 

0.96 

49.93 

0.010 

55 

1000 

1.02 

1.13 

49.92 

0.011 

Table  27: 

Performance  of  MLE  for  a  =  5,  v0  = 

20,  r  =  50,  and  A  =  10,000 

K 

%  Increase 

a&l 

RRlSCzv) 

Tv 

Std.error( rv ) 

24 

20 

35.10 

27.55 

64.27 

0.236 

28 

40 

15.27 

10.44 

53.07 

0.100 

32 

60 

8.49 

6.77 

50.78 

0.067 

36 

80 

5.35 

4.97 

50.11 

0.050 

40 

100 

3.89 

4.73 

49.77 

0.047 

44 

120 

3.11 

4.29 

49.67 

0.043 

48 

140 

2.53 

4.27 

49.53 

0.042 

52 

160 

2.15 

3.92 

49.53 

0.039 

56 

180 

1.90 

3.66 

49.52 

0.036 

60 

200 

1.72 

3.54 

49.53 

0.035 

80 

300 

1.29 

3.03 

49.61 

0.030 

100 

400 

1.14 

2.21 

49.75 

0.022 

120 

500 

1.08 

1.93 

49.82 

0.019 

140 

600 

1.04 

1.41 

49.88 

0.014 

160 

700 

1.03 

1.30 

49.91 

0.013 

180 

800 

1.02 

0.85 

49.94 

0.009 

200 

900 

1.01 

0.60 

49.96 

0.006 

220 

1000 

1.01 

0.39 

49.98 

0.004 
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Table  28:  Performance  of  MLE  for  a  =  10,  v0=l,  r  =  50,  and  N  =  10,000 


Va 

%  Increase 

A&L 

R  fcte(zv) 

|<  ^ 

Std.error(zv ) 

1.2 

20 

80.75 

86.03 

105.81 

0.655 

1.4 

40 

47.12 

34.74 

68.55 

0.294 

1.6 

60 

29.33 

18.31 

58.14 

0.164 

1.8 

80 

20.04 

12.07 

53.93 

0.114 

2 

100 

14.68 

8.91 

52.31 

0.086 

2.2 

120 

10.94 

7.26 

51.35 

0.071 

2.4 

140 

8.70 

6.38 

50.65 

0.063 

2.6 

160 

7.11 

5.35 

50.47 

0.053 

2.8 

180 

5.76 

4.99 

50.18 

0.050 

3 

200 

4.88 

4.83 

49.96 

0.048 

4 

300 

2.70 

3.84 

49.65 

0.038 

5 

400 

1.91 

3.38 

49.61 

0.034 

6 

500 

1.54 

3.14 

49.62 

0.031 

7 

600 

1.33 

3.02 

49.60 

0.030 

8 

700 

1.22 

2.42 

49.71 

0.024 

9 

800 

1.15 

2.16 

49.77 

0.022 

10 

900 

1.10 

1.98 

49.79 

0.020 

11 

1000 

1.07 

1.61 

49.84 

0.016 

Table  29: 

Performance  of  MLE  for  a  =  10,  v0 

=  5,  r  =  50,  and  A  =  10,000 

R, 

%  Increase 

a&l 

R&S(z>) 

i<  ^ 

Std.error(  zv ) 

6 

20 

51.07 

35.04 

68.49 

0.298 

7 

40 

22.41 

11.77 

53.68 

0.112 

8 

60 

11.93 

7.02 

50.98 

0.070 

9 

80 

7.28 

5.22 

50.26 

0.052 

10 

100 

4.91 

4.30 

49.89 

0.043 

11 

120 

3.66 

4.02 

49.69 

0.040 

12 

140 

2.86 

3.80 

49.64 

0.038 

13 

160 

2.36 

3.78 

49.53 

0.038 

14 

180 

2.07 

3.35 

49.60 

0.033 

15 

200 

1.79 

3.37 

49.54 

0.033 

20 

300 

1.26 

2.71 

49.65 

0.027 

25 

400 

1.11 

2.08 

49.78 

0.021 

30 

500 

1.04 

1.30 

49.89 

0.013 

35 

600 

1.02 

1.09 

49.92 

0.011 

40 

700 

1.01 

0.79 

49.95 

0.008 

45 

800 

1.00 

0.51 

49.97 

0.005 

50 

900 

1.00 

0.47 

49.98 

0.005 

55 

1000 

1.00 

0.70 

49.98 

0.007 
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Table  30:  Performance  of  MLE  for  a  =  10,  v0=20,  r  =  50,  and  N  =  10,000 


K 

%  Increase 

A&L 

RRlS(rv) 

'i  >1 

Std.error( rv ) 

24 

20 

34.94 

20.11 

58.28 

0.183 

28 

40 

12.80 

7.68 

51.05 

0.076 

32 

60 

6.41 

5.33 

49.94 

0.053 

36 

80 

3.97 

4.24 

49.67 

0.042 

40 

100 

2.75 

3.72 

49.58 

0.037 

44 

120 

2.13 

3.39 

49.57 

0.034 

48 

140 

1.77 

3.11 

49.59 

0.031 

52 

160 

1.54 

3.14 

49.57 

0.031 

56 

180 

1.38 

2.99 

49.60 

0.030 

60 

200 

1.28 

2.74 

49.64 

0.027 

80 

300 

1.07 

1.75 

49.80 

0.017 

100 

400 

1.02 

1.32 

49.90 

0.013 

120 

500 

1.01 

0.75 

49.95 

0.008 

140 

600 

1.00 

0.35 

49.98 

0.004 

160 

700 

1.00 

0.29 

49.99 

0.003 

180 

800 

1.00 

0.18 

49.99 

0.002 

200 

900 

1.00 

0.06 

50.00 

0.001 

220 

1000 

1.00 

0.04 

50.00 

0.000 

Table  31: 

Performance  of  MLE  for  a  =  50,  v0  =  1, 

r  =  50,  and  A  =  10,000 

K 

%  Increase 

a&l 

R  RtS(zv) 

Tv 

Std.error(  zv ) 

1.2 

20 

117.43 

103.96 

115.84 

0.805 

1.4 

40 

64.79 

35.47 

67.90 

0.306 

1.6 

60 

39.11 

18.26 

57.18 

0.168 

1.8 

80 

25.67 

11.56 

53.28 

0.111 

2 

100 

17.80 

8.47 

51.84 

0.083 

2.2 

120 

13.02 

6.53 

51.05 

0.064 

2.4 

140 

9.96 

5.52 

50.57 

0.055 

2.6 

160 

7.81 

5.04 

50.29 

0.050 

2.8 

180 

6.35 

4.68 

50.06 

0.047 

3 

200 

5.26 

4.71 

49.87 

0.047 

4 

300 

2.68 

3.53 

49.65 

0.035 

5 

400 

1.81 

3.14 

49.60 

0.031 

6 

500 

1.44 

2.93 

49.65 

0.029 

7 

600 

1.23 

2.61 

49.67 

0.026 

8 

700 

1.13 

1.87 

49.78 

0.019 

9 

800 

1.07 

1.60 

49.83 

0.016 

10 

900 

1.04 

1.53 

49.86 

0.015 

11 

1000 

1.02 

1.29 

49.89 

0.013 
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Table  32:  Performance  of  MLE  for  a  =  50,  v0=5,  r  =  50,  and  TV  =  10,000 


K 

%  Increase 

A&L 

R&S(z>) 

|<  Is* 

Std.error(  rv ) 

6 

20 

75.77 

30.75 

63.39 

0.277 

7 

40 

27.62 

9.54 

51.92 

0.094 

8 

60 

12.84 

5.63 

50.42 

0.056 

9 

80 

7.12 

4.27 

49.92 

0.043 

10 

100 

4.45 

3.54 

49.77 

0.035 

11 

120 

3.14 

3.36 

49.65 

0.033 

12 

140 

2.38 

3.26 

49.60 

0.032 

13 

160 

1.91 

3.11 

49.59 

0.031 

14 

180 

1.63 

3.01 

49.59 

0.030 

15 

200 

1.41 

2.67 

49.64 

0.027 

20 

300 

1.08 

1.61 

49.82 

0.016 

25 

400 

1.01 

0.95 

49.93 

0.010 

30 

500 

1.00 

0.33 

49.98 

0.003 

35 

600 

1.00 

0.08 

50.00 

0.001 

40 

700 

1.00 

0.05 

50.00 

0.001 

45 

800 

1.00 

0.02 

50.00 

0.000 

50 

900 

1.00 

0.00 

50.00 

0.000 

55 

1000 

1.00 

0.01 

50.00 

0.000 

Table  33: 

Performance  of  MLE  for  a 

II 

LT\ 

o 

II 

20,  r  =  50,  and  A  =  10,000 

%  Increase 

A&L 

Rfcte(Tv) 

^  >i 

Std.error(  zv ) 

24 

20 

28.56 

10.87 

52.56 

0.106 

28 

40 

7.76 

4.54 

49.94 

0.045 

32 

60 

3.33 

3.50 

49.61 

0.035 

36 

80 

1.96 

3.13 

49.56 

0.031 

40 

100 

1.45 

2.97 

49.60 

0.029 

44 

120 

1.21 

2.57 

49.66 

0.026 

48 

140 

1.10 

2.37 

49.73 

0.024 

52 

160 

1.04 

1.63 

49.84 

0.016 

56 

180 

1.02 

1.10 

49.91 

0.011 

60 

200 

1.01 

0.76 

49.94 

0.008 

80 

300 

1.00 

0.12 

49.99 

0.001 

100 

400 

1.00 

0.05 

50.00 

0.001 

120 

500 

1.00 

0.00 

50.00 

0.000 

140 

600 

1.00 

0.00 

50.00 

0.000 

160 

700 

1.00 

0.00 

50.00 

0.000 

180 

800 

1.00 

0.00 

50.00 

0.000 

200 

900 

1.00 

0.00 

50.00 

0.000 

220 

1000 

1.00 

0.00 

50.00 

0.000 
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Appendix  B:  Cardinality  and  Coverage  Surface  Plots  for  Step  Change  in  Mean 


Figure  17:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  1,  v0=l,  r=50,  and  N  =  10,000. 


Figure  18:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  1,  v0=5,  r=50,  and  N  =  10,000. 
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Cardinality 

Figure  19:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  1,  v0=20,  r=50,  and  N  =  10,000. 


Figure  20:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  5,  v0=l,  r=50,  and iV  =  10,000. 
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Figure  21:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  5,  v0  =5,  r=50,  and  N  =  10,000. 


Cardinality 


Figure  22:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  5,  vo=20,  r=50,  and  iV  =  10,000. 
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Figure  23:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a-  10,  v0=l,  r=50,  and  N  =  10,000. 


Figure  24:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  vQ,  and 

Reference  Value  D.  a  =  10,  v0  =5,  r=50,  and  N  =  10,000. 
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Figure  25:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a-  10,  v0=20,  r=50,  and  N  =  10,000. 


Figure  26:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D.  a  =  50,  v0=l,  r=50,  and N  =  10,000. 
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Figure  27:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D .  a-  50,  v0=5,  r=50,  and  N  =  10,000. 
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Figure  28:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Set  Cardinality,  Coverage,  Percent  Increase  from  v0,  and 

Reference  Value  D .  a  =  50,  v0=20,  r=50,  and  TV  =  10,000. 
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Appendix  C:  Simulation  Data  for  Step  Change  in  Overdispersion  Parameter 


Table  34:  Performance  of  MLE  for  v  =  5,  a0=  1,  r=50,  and  A  =  10,000 


«« 

%  Decrease 

aSl 

R&S(rff) 

$ 

|<  is* 

Std.error(  r« ) 

0.9 

10 

181.60 

235.80 

210.09 

1.731 

0.8 

20 

97.94 

117.52 

127.45 

0.884 

0.7 

30 

59.05 

62.93 

89.00 

0.494 

0.6 

40 

38.45 

36.12 

69.73 

0.303 

0.5 

50 

26.64 

22.19 

59.73 

0.200 

0.4 

60 

18.84 

14.33 

53.83 

0.138 

0.3 

70 

13.77 

9.95 

50.64 

0.099 

0.2 

80 

10.04 

7.05 

49.54 

0.070 

0.1 

90 

7.38 

5.97 

48.63 

0.058 

Table  35: 

Performance  of  MLE  for  v  =  5,  a0  = 

5,  r=50,  and  N  =  10,000 

aa 

%  Decrease 

A&L 

RSlS  (Ta) 

Ta 

Std.error(  r« ) 

4.5 

10 

218.20 

290.90 

248.63 

2.125 

4 

20 

140.78 

179.69 

171.56 

1.323 

3.5 

30 

89.11 

104.73 

120.50 

0.774 

3 

40 

60.38 

65.26 

91.97 

0.500 

2.5 

50 

40.36 

39.97 

72.80 

0.328 

2 

60 

27.86 

25.10 

62.02 

0.220 

1.5 

70 

18.77 

15.03 

55.04 

0.142 

1 

80 

12.22 

9.54 

51.54 

0.094 

0.5 

90 

7.36 

6.21 

49.78 

0.062 

Table  36: 

Performance  of  MLE  for  v  =  5,  a0  = 

10,  r=50,  and  N  =  10,000 

aa 

%  Decrease 

a&l 

R&S(r  a) 

ta 

Std.error(  ra ) 

9 

10 

246.97 

327.69 

276.17 

2.371 

8 

20 

175.21 

226.82 

205.20 

1.654 

7 

30 

124.48 

153.54 

155.15 

1.119 

6 

40 

85.02 

98.24 

115.80 

0.730 

5 

50 

59.27 

63.04 

90.23 

0.485 

4 

60 

40.55 

39.48 

72.81 

0.322 

3 

70 

26.70 

23.51 

60.72 

0.209 

2 

80 

16.92 

13.65 

53.77 

0.131 

1 

90 

9.26 

7.75 

50.18 

0.078 
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Table  37:  Performance  of  MLE  for  v  =  20,  a0  =  1,  r=50,  and  N  =  10,000 


aa 

%  Decrease 

A&L 

R^S(ra) 

$ 

|<  Is* 

Std.error(  r« ) 

0.9 

10 

151.83 

198.58 

184.50 

1.461 

0.8 

20 

76.51 

89.16 

109.05 

0.668 

0.7 

30 

43.77 

46.07 

77.89 

0.367 

0.6 

40 

27.73 

25.38 

63.51 

0.215 

0.5 

50 

18.87 

15.70 

56.31 

0.144 

0.4 

60 

13.30 

10.27 

52.61 

0.099 

0.3 

70 

9.46 

7.22 

50.53 

0.072 

0.2 

80 

6.77 

5.68 

49.51 

0.057 

0.1 

90 

4.76 

4.57 

48.98 

0.045 

Table  38: 

Performance  of  MLE  for  v  =  20,  a0 

=  5,  r=50,  and  N  =  10,000 

%  Decrease 

aSl 

R®lS(Ta) 

l<  ^ 

Std.error(  r« ) 

4.5 

10 

177.80 

235.03 

209.76 

1.724 

4 

20 

93.53 

114.08 

126.01 

0.851 

3.5 

30 

55.23 

61.16 

88.53 

0.475 

3 

40 

35.67 

36.15 

70.57 

0.297 

2.5 

50 

23.46 

21.45 

60.60 

0.186 

2 

60 

15.53 

13.51 

54.73 

0.127 

1.5 

70 

10.52 

9.17 

52.00 

0.090 

1 

80 

6.77 

6.23 

50.51 

0.062 

0.5 

90 

3.99 

4.11 

49.86 

0.041 

Table  39: 

Performance  of  MLE  for  v  =  20,  a0 

=  10,  r=50,  and  A  =  10,0C 

au 

%  Decrease 

aSl 

Rfe(r„) 

l<  ^ 

Std.error(  r« ) 

9 

10 

192.08 

252.41 

223.78 

1.831 

8 

20 

112.28 

140.27 

144.19 

1.039 

7 

30 

66.85 

76.22 

99.57 

0.579 

6 

40 

43.35 

45.26 

77.37 

0.360 

5 

50 

28.19 

26.67 

64.21 

0.226 

4 

60 

19.11 

16.79 

57.04 

0.152 

3 

70 

12.62 

10.87 

52.88 

0.105 

2 

80 

8.02 

7.47 

50.93 

0.074 

1 

90 

4.44 

4.76 

50.01 

0.048 
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Appendix  D:  Cardinality  and  Coverage  Surface  Plots  for  Step  Change  in 

Overdispersion  Parameter 


Cardinality 


Figure  29:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=5,  a0= 1,  r=50,  and  A  =  10,000. 


Cardinality 


Figure  30:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=5,  a0= 5,  r=50,  and  N  =  10,000. 
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Figure  31:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=5,  a0= 10,  r=50,  and  N  =  10,000. 


Figure  32:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=20,  a0= 1,  r= 50,  and  =  10,000. 


89 


Cardinality 


Figure  33:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=20,  a0=  5,  r= 50,  and N  =  10,000. 


Figure  34:  Surface  Plot  Obtained  from  Confidence  Set  Estimator  Showing  Estimated 
Relationships  Between  Cardinality,  Coverage,  Percent  Decrease  from  a0 ,  and  Reference 

Value  D.  v=20,  a0  =  10,  r=50,  and  iV  =  10,000. 
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